-2.1cmHessian heat kernel estimates and
Calderón-Zygmund inequalities
on complete Riemannian manifolds
* Corresponding author 2010 Mathematics Subject
Classification. Primary: 35K08; Secondary: 35J10, 35J30, 47G40.
Key words and phrases.  Brownian motion, Heat
semigroup, Heat kernel, Bismut type Hessian representations,
Hessian formula, Calderón-Zygmund inequalities, the
Hardy-Littlewood maximal function, good-λ inequalities,
Kato class.  This work has been project is supported by
the National Natural Science Foundation of China (Grant
No. 12071431). The third author has been supported by the Fonds
National de la Recherche Luxembourg (project GEOMREV

1]Jun Cao 2]Li-Juan Cheng* 3]Anton Thalmaier

[1]Department of Applied Mathematics, Zhejiang
University of Technology,Hangzhou 310023, The People's Republic of China

[2]School of  Mathematics, Hangzhou Normal UniversityHangzhou 311123, The People's Republic of China

[3]Department of Mathematics, University of Luxembourg, Maison du
Nombre,L-4364 Esch-sur-Alzette, Luxembourg

We address some fundamental questions
about geometric analysis on Riemannian manifolds. The
$L^p$-Calderón-Zygmund inequality is one of the cornerstones in
the regularity theory of elliptic equations, and it has been
asked under which geometric conditions it holds for a reasonable class
of non-compact Riemannian manifolds, and to what extent
assumptions on the derivative of curvature and on the
injectivity radius of the manifold are necessary.
In the present paper, for $1<p<2$, we give
a positive answer for the validity of the $L^p$-Calderón-Zygmund
inequality on a Riemannian manifold assuming only a lower
bound on the Ricci curvature. It is well known that  this  alone is
not sufficient for $p>2$.
In this case we complement the study of Güneysu-Pigola (2015) and
derive sufficient geometric criteria for the validity of the
Calderón-Zygmund inequality under additional Kato class bounds
on the Riemann curvature tensor and the covariant derivative of
Ricci curvature. Bounds in the Kato class are
integral conditions and much weaker than pointwise
bounds. Throughout the proofs, probabilistic tools, like Hessian
formulas and Bismut type representations for heat semigroups, play
a significant role.

§ INTRODUCTION

The Hessian operator $\Hess$, which contains all information of the second order derivatives, is a fundamental object in the second order smooth analysis.
To study the Hessian operator, one usually needs
the Calderón-Zygmund inequality that controls
$\Hess$ by the simpler Laplacian operator $\Delta$
(see [6, 20, 35, 26, 36] and the references therein for its wide applications particularly in the regularity theory of elliptic equations).

In the Euclidean space $\mathbb{R}^d$, the Calderón-Zygmund inequality says that
for any $p\in (1,\infty)$ and $u\in C_c^{\infty}(\R^n)$, it holds

\begin{align}\label{eqn-CZes}
\|\Hess u\|_{L^p(\R^d)}\leq C\| \Delta u\|_{L^p(\R^d)},
\end{align}

where $C=C(d,p)>0$ is a constant depending only on $d$ and $p$.
The inequality (<ref>) is first proved by Calderón and Zygmund
via their seminal
theory of singular integral operators based on the explicit
representation of the Green kernel of the Laplacian.
A remarkable consequence of inequality (<ref>) is the fact that,
in the Euclidean space, the Sobolev norms

\begin{align*}
\|u\|_{W^{2,p}(\R^d)}&=\|u\|_{L^p(\R^d)}+\|\nabla u\|_{L^p(\R^d)}+\|\Hess u\|_{L^p(\R^d)},\\
\|u\|_{\widetilde{W}^{2,p}(\R^d)}&=\|u\|_{L^p(\R^d)}+\|\Delta u\|_{L^p(\R^d)}
\end{align*}

are equivalent on $C_c^{\infty}(\R^n)$ (see [36]).

The  Calderón-Zygmund inequality (<ref>)
extends to second order uniformly elliptic operators $L=-\div A \nabla$ with variable coefficients $A$ on domains $\Omega\subset \R^d$ without any boundary conditions (see e.g. [20]). In this setting, we have the following local Calderón-Zygmund inequality that
given any domains $\Omega_1\Subset  \Omega$,  $p \in (1,\infty)$ and
$u\in C_c^{\infty}(\Omega)$, it holds

\begin{align}\label{eqn-CZlv}
\|\Hess  u\|_{L^p(\Omega_1)}\leq C\Big(\|u\|_{L^p(\Omega)}+\|Lu\|_{L^p(\Omega)}\Big),
\end{align}

where $C=C(\Omega_1,\Omega,d,p,A)>0$ is a constant
depending on $\Omega_1, \Omega, p,d$ and the elliptic coefficients of $A$. The proof of (<ref>) is based on (<ref>) and a perturbation argument of $A$. In particular, if $A$ is a constant matrix,
then one can get rid of
the term $\|u\|_{L^p(\Omega)}$ in (<ref>).

A further step was taken by Güneysu and Pigola in [26] where they considered the following global Calderón-Zygmund inequality
on Riemannian manifolds $M$ of the form that
for any $p\in (1,\infty)$ and $u\in C_c^\infty(M)$,
it holds

$$\|\Hess u\|_{L^p(M)}\leq C_1 \|u\|_{L^p(M)}+C_2\|\Delta u\|_{L^p(M)}, \quad \eqno{{\bf CZ}(p)}$$

where $\Delta$ is the Laplace-Beltrami operator on $M$ and $C_1$,
$C_2$ are two positive constants.  It is known that such an inequality
CZ($p$) can hold or fail depending on $p$ and the geometry of
$M$. We give a brief summary of the state of the art on this subject
(see [36] for a more detailed survey, as well as
[27]).  For an example of a Riemannian manifold of positive
sectional curvature where CZ($p$) fails for large values of $p$,
see [33], as well as
[30] for a general result in this direction.

To begin with, we make some conventions on the notation.  Throughout
this paper, let $(M,g)$ be a complete non-compact connected
$d$-dimensional Riemannian manifold, $\nabla$ the Levi-Civita
covariant derivative, $\Hess=\nabla d$ the Hessian operator on
functions, and $\mu$ the Riemannian volume measure on $M$.  We denote
by $|\cdot|$ the norm in the tangent space, and by $\|\cdot\|_p$ the
norm in $L^p(M,\mu)$ for $1\leq p\leq \infty$. The Laplace-Beltrami
operator $\Delta$ acting on functions and forms, is understood as
self-adjoint positive operator on $L^2(\mu)$.

If $p=2$ and there is a bound $\Ric\geq -K$ of the Ricci curvature of
$M$ for some constant $K>0$, then it is well-known that CZ($2$)
is a straightforward consequence of Bochner's identity, see the
Appendix. The extension of CZ($p$) from $p=2$ to an arbitrary
$p\in (1,\infty)$ is much more involved and an intriguing
problem. Inspired by the proof of (<ref>), a possible way to
establish ${\bf CZ}(p)$ is a similar potential theoretical approach that
represents $u$ via the Green kernel of the Laplace-Beltrami
operator $\Delta$ (see [35, 41, 34]). Although
the method works pretty well for many related questions of the
associated Poisson equation, it has the drawback that it applies only
under some restrictive conditions such as compactness of $M$ or nonnegative Ricci

To overcome this drawback, Güneysu and Pigola [26, 36]
introduced two methods of proof for CZ$(p)$ avoiding the use of
Green kernel. The first method is based on a gluing
procedure that connects the local consideration to the global
result. To be precise,  if $M$ has bounded Ricci curvature and a
strictly positive injectivity radius, then it is proved in
[26] that the Calderón-Zygmund
inequality CZ$(p)$ holds for any $p\in (1,\infty)$ with
implicit constants depending on $p$, $\|{\Ric}\|_{L^\infty}$, the
dimension $d$  and the injectivity radius. Moreover, if
$p>\max\{2,~d/2\}$ and  $M$ has bounded sectional curvature, then CZ$(p)$
also holds (see <cit.>). The second method,
called the functional analytic method, uses boundedness results for
the covariant Riesz transform for $1<p<2$ from [40].

Let us sketch the main idea of the second method in [26].
Inequality CZ($p$) is usually reduced to the existence of
positive constants $C$ and $\sigma$ such that

\begin{align*}
\left\|\,|\Hess (\Delta+\sigma)^{-1}u|\,\right\|_p\leq C\,\|u\|_p,
\end{align*}

which is equivalent to

\begin{align}\label{eqn-equiv1}
\left\|\,|\nabla (\Delta_1+\sigma)^{-1/2}
\circ d(\Delta_0+\sigma)^{-1/2}u|\,\right\|_p\leq C\,\|u\|_p.
\end{align}

Here and hereafter, we write $\|\cdot\|_{p}:=\|\cdot\|_{L^p(M)}$ for simplicity.
The problem is thus reduced to the study of conditions for boundedness
of the classical Riesz transform $d(\Delta_0+\sigma)^{-1/2}$ on
functions and boundedness of the covariant Riesz transform
$\nabla (\Delta_1+\sigma)^{-1/2}$ on one-forms.  This approach in
[26] however is restricted to $p\in (1,2)$; the constants
$C_1, C_2$ in CZ($p$) depend on dimension $d$, $p$, $\|R\|_{\infty}$,
$\|\nabla R\|_{\infty}$ and on the constants $D, \delta$ from the
following local volume doubling assumption: there are constants $C>0$,
$0\leq \delta < 2$ such that

\begin{align}\label{eqn-GD}
V(x,tr)\leq Ct^d\e^{Ctr}V(x,r) \tag{\bf LD}
\end{align}

for all $x\in M$, $r>0$ and $t\geq 1$, where $V(x,r):=\mu(B(x,r))$.

Very recently, Baumgarth, Devyver and Güneysu [4]
studied the covariant Riesz transform on $j$-forms. Their results can
be applied to CZ($p$) when $1<p<2$ requiring the same curvature
conditions as in [26] but without the local volume doubling
assumption made in [26]. This comes from the fact that
(<ref>) already holds if $\Ric\geq -K$ for some $K\geq 0$, as
can be seen by the Bishop-Gromov comparison theorem and the well-known
formula for the volume of balls in hyperbolic space.  However, it
seems difficult to establish CZ($p$) for $p>2$ in this way, as
when trying to extend the machinery of [3] to
$L^p$-boundedness of covariant Riesz transform for $p>2$, the local
Poincaré inequality, used explicitly in [3], does not
make sense on differential forms.

The observations above raise the following questions:

* In the case $1<p<2$, is it possible to weaken the
assumptions on the Riemann curvature tensor $\|R\|_{\infty}$ and
$\|\nabla R\|_{\infty}$, e.g., replacing them by Ricci curvature

* Without a lower control on the injectivity radius, under
which conditions on the manifold $M$, the inequality CZ($p$)
holds on $M$ for $p>2$?

In order to answer the above two questions affirmatively, we develop a
new functional analytic method of proof for CZ$(p)$ that works
for all $p\in (1,\infty)$. Unlike the second method used in
[26, 36] that reduces CZ($p$) to the boundedness
(<ref>), our method makes use of the observation that
CZ($p$) is equivalent to the $L^p$-boundedness of the following

\begin{align*}
\Hess(\Delta+\sigma)^{-1}=\int_0^{\infty}\e^{-\sigma t}\Hess P_t\,dt,
\end{align*}

where $P_t$ denotes the heat semigroup generated by $-\Delta$. Based on this observation, our strategy
to solve the above two questions is first establishing some Hessian heat kernel estimates and then bridging from Hessian heat kernel estimates to CZ($p$).

However, references on Hessian heat kernel estimates are sparsely to find in the literature,
especially when compared to the situation of heat kernel estimates and
gradient estimates for the heat kernel. Therefore, in the following section, we invest some
effort in deriving sharp Hessian heat kernel estimates first.
It turns out that there is a big difference
between the treatment of the cases $p<2$ and $p>2$.
In particular, if $p\in (1,2)$, we only assume that $M$ satisfies the following curvature condition:

\begin{align}\label{eqn-Ric}
\Ric\geq -K    \tag{\bf Ric}
\end{align}

for some $K\geq 0$. The curvature condition (<ref>) then give rise to the local Gaussian upper bounds of the heat kernel and the local volume doubling condition (LD) which enable us to deduce a series of
$L^2$ weighted off-diagonal estimates of the Hessian heat kernel (see Section <ref> below). These estimates
together with the classical argument of Calderón-Zygmund decomposition (see [12]) allow to derive our first main result.

Let $(M,g)$ be a complete Riemannian manifold satisfying
(<ref>). Let $1<p<2$ be fixed. Then there exists a constant
$\sigma>0$ such that the operator $\Hess (\Delta+\sigma)^{-1}$ is
bounded in $L^p$, i.e. CZ($p$) holds.

Theorem <ref> answers question $1$ affirmatively. Comparing Theorem <ref> with existing results on CZ($p$),
it should be pointed out first that the result is valid without any
injectivity radius assumptions and secondly, instead of boundedness of
$\|R\|_{\infty}$ and $\|\nabla R\|_{\infty}$, only a lower bound of
$\Ric$ is needed.  In view of the argument for $p=2$, the assumption (<ref>) seems close to

For the case $p>2$,  it is well known
[36, 33, 30] that the curvature condition (<ref>)
alone is not enough for CZ($p$) to hold.
Hence more geometric information about the manifold is
required in this case for the validity of CZ($p$).

To formulate appropriate geometric conditions, we introduce some  probabilistic quantities.
Denote by $X_t$ the diffusion process generated by $-\Delta$, which is
assumed to be non-explosive ($M$ is stochastically complete). Then,
for $f\in \mathcal{B}_b(M)$,

\begin{align*}
\E^x[f(X_t)]=P_tf(x)=(\e^{-\Delta t}f)(x).
\end{align*}

A Borel function $K\colon M\rightarrow \R$ is said to be in the Kato
class $\mathcal{K}(M)$ of $M$, if

\begin{align}\label{eqn-Katoc}
\lim_{t\rightarrow 0+}\sup_{x\in M}\int_0^t \E^x\big[|K(X_s)|\big]\,ds=0. \tag{{\bf K}}
\end{align}

Obviously, $\mathcal{K}(M)$ is a linear space and
$\mathcal{K}(M)\subset L_{\rm loc}^1(M)$.
The Kato class has been introduced in [31] on Euclidean space and then applied
to investigate singular potentials, see for instance
[1, 38, 22, 37].
Concerning criteria for functions to be in the Kato class the reader
may consult [24, 23, 28].  Note that for dimension $d\geq 2$, if (<ref>) holds, then the heat kernel has a local on-diagonal estimate which implies $L^p(M)+L^{\infty}(M)\subset \mathcal{K}(M)$ for $p>d/2$ by <cit.>.

To deal with the case $p>2$,
our conditions are
given in terms of Kato bounds on the  geometric quantities.

Condition (H)   $\Ric\geq -K$ for some $K\geq 0$
and there exist $K_1,K_2\in \mathcal{K}(M)$ such that

\begin{align}\label{eqn-H}
|R|^2(x)\leq K_1(x)\quad\text{and}\quad  \big|\nabla \Ric^{\sharp}+d^*R\big|^2(x)\leq K_2(x),\quad
x\in M, \tag{\bf H}
\end{align}

where  for $x\in M$ and $ v_1,v_2,v_3\in T_xM$,

\begin{align*}
&|R|(x)=\sup \Big\{|R^{\#,\#}(v_1,v_2)|_\HS(x)\colon\ v_1,v_2\in T_xM,\ |v_1|\leq 1,  \ |v_2|\leq 1 \Big\},
\end{align*}

\begin{align*}
& \langle d^*R(v_1,v_2), v_3 \rangle=\langle (\nabla_{v_3}\Ric^{\sharp})(v_1), v_2 \rangle-
\langle (\nabla_{v_2}\Ric^{\sharp})(v_3), v_1 \rangle,
\end{align*}

with $R^{\#,\#}(v_1,v_2)=R(\cdot,v_1,v_2,\cdot)$, the  curvature tensor $R$ and $\Ric^{\sharp}(v)=\Ric(\cdot, v)^\sharp$ for $v\in T_xM$.

Under the condition ${\bf (H)}$, we are able to prove the
following key pointwise inequalities for the Hessian of the semigroup that for any $t>0$,
$f\in C_c^\infty(M)$ and $x\in M$, it holds

\begin{align}\label{eqn-pihe}
|\Hess P_tf|(x)\leq \e^{2Kt}P_t|\Hess f|(x) + C\e^{(2K+\theta)t}(P_{t}|\nabla f|^2)^{1/2}(x),
\end{align}

\begin{align}\label{eqn-Bismut-He}
&t\,|\Hess P_tf|(x)\leq C(1+\sqrt{t})\,\e^{(2K+\theta)t}\lf(P_t|f|^2\r)^{1/2}(x)
\end{align}

for some constants $C,\theta>0$ (see Propositions <ref> and
<ref> below).  Both of them are proved by using some
probabilistic tools established in
[18, 19, 32, 39, 2]. More precisely, inequality
(<ref>) is proved by a stochastic approach based on a second
order derivative formula for the heat semigroup (see
(<ref>) below); inequality (<ref>) is
established by means of Bismut-type representation formulas for the
Hessian of heat semigroups which were first proved by Elworthy and Li
(see [18, 19]).

The inequalities (<ref>) and (<ref>) enable us to
establish a series of pointwise Hessian heat kernel estimates.
Based on these results, we can give an affirmative answer to the question $2$ above.

Let $M$ be a complete Riemannian manifold satisfying
(<ref>). Let $p>2$ be fixed. Then there exists a constant
$\sigma>0$ such that the operator $\Hess\,(\Delta+\sigma)^{-1}$ is
bounded in $L^p$, i.e. CZ$(p)$ holds.

Theorem <ref> gives in particular an answer to the open
question in [26] about sufficient conditions for CZ($p$)
when $p>2$ in the absence of control of the injectivity radius. It is
worth mentioning that compared with the sufficient conditions even in
the case $1<p<2$ (<cit.>,
<cit.>), in Theorem <ref> the
geometric quantities $|R|$ and $|\nabla \Ric|$ do not need to be
uniformly bounded on $M$; it is sufficient to have bounds in the Kato
class which is a kind of integral condition.
It should be mentioned that for $p>d/2$ in conjunction with
$p\geq2$ the validity of CZ$(p)$ has been obtained under
strong curvature assumption, namely uniform boundedness of the
sectional curvature tensor, but without conditions involving the
derivative of curvature, see Theorem 5.18 in the survey of
Pigola [36].

As indicated above, the pointwise inequalities (<ref>) and
(<ref>) play a key role in the proof of Theorem
<ref>, which also reflects the difference between our
approach and that of classical local Riesz transform (see
[13, 14, 3]).  In the latter case, the following
domination property that for some positive constants $c_1,c_2$ and $C$

\begin{align}\label{Pt-one}
|\nabla P_tf|\leq C\e^{c_1t}P_{c_2t}|\nabla f|
\end{align}

is indispensable for boundedness of the
Riesz transform in case $p>2$ (see [13, 14]).
Unfortunately, the machinery of [13] cannot work well in our situation since a
Hessian estimate of the type

\begin{align*}
|\Hess P_tf|\leq C\e^{c_1t}P_{c_2t}|\Hess f|
\end{align*}

would be required which however only holds in very specific cases (like flat
manifolds). On the other hand, if we aim at using the techniques from
[3] directly, the main difficulty to deal with is that there is no
suitable Hessian replacement of the local Poincaré inequality which is
heavily used throughout their proof.

To overcome these two obstacles, we take advantage of the techniques
from [3], in particular, the sharp maximal function and
good-$\lambda$ inequalities. By means of these tools and inequality
(<ref>), we observe that if the additional term involving
$(P_t|\nabla f|^2)^{1/2}$ is treated and controlled as an error term, then
boundedness of $\Hess (\Delta+\sigma)$ in $L^p$ can be
established also for $p>2$ with the help of pointwise Hessian estimates
of the heat kernel.  In conclusion, the crucial observation is that the pointwise inequality
(<ref>) may serve as a Hessian replacement of
(<ref>) circumventing the non-availability
of the local Poincaré inequality.

The rest of the paper is organized as follows. In Section 2, we give
various forms of estimates on the Hessian of the heat kernel and on
the corresponding semigroups, which are used throughout the proof of
CZ($p$). In Sections 3 finally, we present proofs for Theorem
<ref> and <ref> respectively. A suitable version of
the localization techniques of [17] is included.

Acknowledgements The authors are indebted to Batu
Güneysu, Stefano Pigola and Giona Veronelli for very helpful
comments on the first draft of this paper.

§ HESSIAN HEAT KERNEL  ESTIMATES

This section is divided into two parts: the first part is on
$L^2$-estimates of the Hessian of heat kernels under the assumption
of a lower Ricci curvature bound; the second part is on the estimates
derived under condition (H) from the Hessian formula and Bismut
type formulas for the Hessian of semigroups, which are established by
means of the techniques from stochastic analysis.

§.§ $L^2$-estimates for the Hessian of heat kernel

In this section, we assume Ricci curvature to be bounded below,
i.e. validity of condition (Ric).  Then, in particular, the
local doubling assumption (<ref>) with respect to $\mu$ holds.
Then for $x$ and $y\in M$, we obviously have
$B(y,\!\sqrt{t})\subset B(x,\sqrt{t}+\rho(x,y))$. Thus if
(<ref>) holds, then

\begin{align}\label{eqn-gdxy}
V(y,\!\sqrt{t})\leq V\Big(x,\sqrt{t}+\rho(x,y)\Big)
\leq C\Big(1+\frac{\rho(x,y)}{\sqrt{t}}\Big)^{d}\exp\left(C(\sqrt{t}
\end{align}

It is well-known that, under a lower Ricci curvature bound, the heat
kernel allows an off-diagonal estimate [15].  The following
lemma gives a pointwise off-diagonal estimate for the heat kernel and
its time derivative.

Assume that (<ref>) holds.  Then for any
$\az\in (0,\frac{1}{4})$, there exist constants $C$ and $C_1>0$
depending on the dimension $d$ and $\alpha$ such that for all
$x,y\in M$ and $t>0$,

\begin{align}\label{heat-kernel-upper-bound}
p_t(x,y)+\left|\frac{\partial p_t}{\partial t}(x,y)\right| \leq \frac{C}{V(y,\!\sqrt{t})}
\exp\left(-\alpha \frac{\rho^2(x,y)}{t}+C_1 Kt\right).
\end{align}

The estimate of $p_t(x,y)$ is an easy consequence of
(<ref>) and <cit.> or <cit.>, where it is proved that for all $x,y\in M$ and

\begin{align}\label{heat-kernel-upper-bound-2}
p_t(x,y)\leq \frac{C}{V(y,\!\sqrt{t})}
\exp\lf(-\alpha \frac{\rho^2(x,y)}{t}+C_1 Kt\r).
\end{align}

The estimate for $|\frac{\partial p_t}{\partial t}(x,y)|$ follows
from that of $p_t(x,y)$ and the analytic property of the semigroup
(see <cit.> or <cit.>).

The following lemma gives weighted $L^2$-integral estimates for the
heat kernel, its gradient and its Laplacian.

Assume that (<ref>) holds.  Let $\az\in (0,\frac{1}{4})$ be
as in Lemma .  For all
$\gamma \in (0,2\alpha)$, $s>0$ and $y\in M$,

\begin{align*}
\int_M \lf[|p_s(x,y)|^2+s|\nabla _x p_s(x,y)|^2+s^2|\Delta_x p_s(x,y)|^2\r]
\e^{\gamma\, \frac{\rho^2(x,y)}{s}}\,\mu(dx)
\leq \frac{C_{\gamma}}{V(y,\!\sqrt{s})}\,\e^{2\,C's},
\end{align*}

where $C_{\gamma}>0$ depends on $\gamma$ and $C'>0$ on $\az$ and

By (<ref>), it is easy to see that for all $\gamma>0$,
$s,t>0$ and $y\in M$, there exist two positive constants
$C_{\gamma}$ (depending on $\gamma$ and the constants in
(<ref>)) and $\tilde{C}$ such that

\begin{align}\label{eqn-ele1}
\int_{\rho(x,y)\geq \sqrt{t}}\e^{-2\gamma\frac{\rho^2(x,y)}{s}}\,\mu(d x)
&\leq \e^{-\gamma t/s}\int_M \e^{-\gamma \frac{\rho^2(x,y)}{s}}\,\mu(dx) \notag\\
&\leq  \e^{-\gamma t/s}\sum_{i=0}^{\infty}V(y,(i+1)\sqrt{s})\e^{-\gamma i^2}\notag\\
&\leq C  \e^{-\gamma t/s}V(y,\!\sqrt{s}) \sum_{i=0}^{\infty}
(i+1)^d \e^{-\gamma i^2}\e^{C(i+1)\sqrt{s}} \notag\\
&\leq  C  \e^{-\gamma t/s}\e^{C\sqrt{s}}V(y,\!\sqrt{s}) \sum_{i=0}^{\infty}
(i+1)^d \e^{-\gamma i^2}\e^{\gamma i^2/2+C^2s/2}\notag\\
& \leq  C  \e^{-\gamma t/s}\e^{C\sqrt{s}+C^2s/2}V(y,\!\sqrt{s}) \sum_{i=0}^{\infty}
(i+1)^d \e^{-\gamma i^2/2}\notag\\
& \leq C_{\gamma} V(y,\!\sqrt{s}) \,\e^{-\gamma t/s} \e^{\tilde{C}s},
\end{align}

where the second inequality comes from condition (<ref>). The
remainder of the proof follows from <cit.>.

We now turn to the estimates for the Hessian of heat kernel. The
following lemma shows that the Hessian of heat semigroup also
satisfies an $L^2$-Gaffney off-diagonal estimate. Note that in the
following discussion the constant $C$ will be different in different
lines without confusion.

Assume that (<ref>) holds.  There exist constants $C,C_2>0$
such that for all $t\in (0,\fz)$, all Borel subsets $E,F\subset M$
with compact closure, and all $f\in L^2(M)$ with supp $f\subset E$,

\begin{align*}
\big\|\mathbbm{1}_Ft\, |\Hess P_tf|\big\|_2\leq C(1+\sqrt{t})\,\exp\left(-\frac{C_2\,\rho^2(E,F)}{t}\right)\|f\|_2.
\end{align*}

Recall the following $L^2$-Gaffney off-diagonal estimate on
one-forms [4]. If (<ref>) holds, then for
$\alpha \in \Gamma_{L^2}(T^*M)$ with support
${\rm supp}(\alpha)\subset E$ and any $s\in (0,1)$, it holds

\begin{align*}
\big\|\mathbbm{1}_F\!\sqrt{s}\,|\nabla \e^{-s\Delta^{(1)}}\alpha|\big\|_2
&\leq C\,(1+\sqrt{s})\,\exp\left(-\frac{c_1\,\rho^2(E,F)}{s}\right)\|\mathbbm{1}_E \alpha\|_2\\
&\leq 2C\exp\left(-\frac{c_1\rho^2(E,F)}{s}\right)\|\mathbbm{1}_E \alpha\|_2
\end{align*}

for some positive constants $C$ and $c_1$.  On the other hand, we
have Gaffney's off-diagonal estimate for $|\nabla P_{t_2}f|$ (see
<cit.>), i.e. for all $f \in L^2(M)$ with support
in $E$ and $u>0$,

\begin{align*}
\big\|\mathbbm{1}_{F}\!\sqrt{u}\,|\nabla P_{u}f|\big\|_2\leq C\exp\left(-\frac{c_2\rho^2(E,F)}{u}\right)\|\mathbbm{1}_E f\|_2
\end{align*}

for some positive constants $C$ and $c_2$.  For $t>0$, denoting by
$\Delta^{(1)}$ the Laplacian on one-forms, we may write

$$\Hess P_tf=\frac{\sqrt{\left(t-(\frac{t}{2}\wedge 1)\right)}\,\nabla \e^{-\big(t-(\frac{t}{2}\wedge 1)\big)\Delta^{(1)}}\left(\sqrt{\frac{t}{2}\wedge 1}\,d P_{\frac{t}{2}\wedge 1}f\right)}{\sqrt{\left(t-(\frac{t}{2}\wedge 1)\right)\big(\frac{t}{2}\wedge 1\big)}}  $$

so that using the composition rule of Gaffney's off-diagonal estimate
(see <cit.>), we obtain

\begin{align*}
\big\|\mathbbm{1}_{E}t\, |\Hess P_{t}f|\big\|_2\leq C\frac{t}{\sqrt{\left(t-(\frac{t}{2}\wedge 1)\right)(\frac{t}{2}\wedge 1)}}\exp\left(-\frac{c\,\rho^2(E,F)}{\max\left\{\frac{t}{2}\wedge 1, t-(\frac{t}{2}\wedge 1)\right\}}\right)\|\mathbbm{1}_E f\|_2,
\end{align*}

for some positive constants $C$ and $c$.  Note that for $0<t/2<1$,

\begin{align*}
\max\left\{\frac{t}{2}\wedge 1, t-\big(\frac{t}{2}\wedge 1\big)\right\}=\max\left\{\frac{t}{2}, t-\frac{t}{2}\right\}=\frac{t}{2};\quad
\frac{t}{\sqrt{\left(t-\big(\frac{t}{2}\wedge 1\big)\right)(\frac{t}{2}\wedge 1)}}=2,
\end{align*}

and that for $t/2\geq 1$,

\begin{align*}
\max\left\{\frac{t}{2}\wedge 1, t-\big(\frac{t}{2}\wedge 1\big)\right\}=\max\big\{1, t-1\big\}=t-1\leq t;\quad
\frac{t}{\sqrt{\left(t-\big(\frac{t}{2}\wedge 1\big)\right)(\frac{t}{2}\wedge 1)}}=\frac{t}{\sqrt{t-1}}\leq \sqrt{2t}.
\end{align*}

Therefore, we conclude that there exist positive constants $C$ and
$C_2$ such that

\begin{equation*}
\big\|\mathbbm{1}_{F}t\, |\Hess P_{t}f|\big\|_2\leq C\,(1+\sqrt{t})\,\exp\left(-\frac{C_2\rho^2(E,F)}{t}\right) \|\mathbbm{1}_E f\|_2.\qedhere
\end{equation*}

The following proposition gives $L^2$-weighted estimates for the
Hessian of the heat kernel.

Assume that (<ref>) holds.  Fix $\az\in (0,\frac{1}{4})$ as
in Lemma .  Then for all
$\gamma\in (0,2\alpha)$, $s>0$ and $y\in M$, there exists a constant
$C>0$ such that

\begin{align*}
\int_M |\Hess_xp_s(x,y)|^2\,\exp\left(\gamma\,\frac{\rho^2(x,y)}{s}\right)\,\mu(dx)\leq \frac{C(1+Ks)\,\e^{2C's}}{s^2 V(y,\!\sqrt{s})}
\end{align*}

where $C'>0$ is the same constant as in Lemma <ref>.

We begin by integrating Bochner's identity (<ref>) to obtain

\begin{align*}
\frac{1}{2}&\int |\nabla p_s|^2 \,\Delta \e^{\gamma\rho^2(x,y)/s}\,  \mu(dx)
=\frac{1}{2}\int_M (\Delta |\nabla p_s|^2)\,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)\\
&=\int_M \left(-|\Hess p_s|^2_\HS+g(\nabla \Delta p_s, \nabla p_s )-\Ric(\nabla p_s,\nabla p_s)\right)\e^{\gamma\rho^2(x,y)/s}\,\mu(dx),
\end{align*}

which then implies

\begin{align}\label{eq-Hess-upper-bound}
&\int_M |\Hess p_s|_\HS^2\,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx) \\
&=-\frac{1}{2}\int_M |\nabla p_s|^2\Delta \e^{\gamma\rho^2(x,y)/s}\,\mu(dx)
+\int_Mg( \nabla \Delta p_s, \nabla p_s) \,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)\notag\\
\Ric (\nabla p_s, \nabla p_s)\,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)\notag\\
&=-\frac{\gamma}{s}\int_M |\nabla p_s|^2\div\left(\e^{\gamma\rho^2(x,y)/s}\rho(x,y)
{\nabla}\rho(\cdot, y)(x)\right)\,\mu(dx) + \int \Delta p_s\div (\e^{\gamma\rho^2(x,y)/s}\nabla p_s)\,\mu(dx)\notag\\
&\quad-\int \Ric(\nabla p_s, \nabla p_s)\,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)
=:\mathrm{I}_1+\mathrm{I}_2+\mathrm{I}_3 \notag
\end{align}

(with the sign convention $\Delta=\div\nabla$ for the
divergence). Now for any $\beta\in (0,2\az)$, let

\begin{align*}
E(s,y,\beta):=\int_M \lf[|p_s(x,y)|^2+s|\nabla _x p_s(x,y)|^2+s^2|\Delta_x p_s(x,y)|^2|\r]
\e^{\beta\rho^2(x,y)/s}\,\mu(dx).
\end{align*}

Then by condition (<ref>) and Lemma <ref>, we have

\begin{align}\label{eqn-esI3}
\mathrm{I}_3&\le K\int_M |\nabla p_s|^2\,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)
\le \frac{K}{s} E(s,y,\gz)\le \frac K s \frac{C_{\gamma}}{V(y,\!\sqrt{s})}\,\e^{2 C's}.
\end{align}

For $\mathrm{I}_2$, using the fact $|\nabla\rho|\le 1$, we deduce
from Lemma <ref> again that

\begin{align}\label{eqn-esI2}
\mathrm{I}_2&=\int_M (\Delta p_s)^2\,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)
-\frac{2\gamma}{s}\int_M (\Delta p_s)\,\e^{\gamma\rho^2(x,y)/s} \rho(x,y)\langle\nabla \rho, \nabla p_s \rangle \,\mu(dx) \notag \\
&\leq 2\int_M (\Delta p_s)^2\,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)+
\frac{\gamma^2}{s^2}\int_M \e^{\gamma\rho^2(x,y)/s}\rho^2(x,y)|\nabla p_s|^2\,\mu(dx)\notag\\
&\le  \frac{C_\gz}{s^2}E(s,y,\gz')\le \frac 1 {s^2} \frac{C_{\gamma'}}{V(y,\!\sqrt{s})}\,\e^{2 C's},
\end{align}

where $\gamma<\gamma'<2\alpha$ with $\gamma'-\gamma$ sufficiently

To bound $\mathrm{I}_1$, we write

\begin{align}\label{eqn-esI1-1}
\frac s{\gamma}\mathrm{I}_1
&=\int_M |\nabla p_s|^2\,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx) +2
\int_M |\nabla p_s|^2\,\e^{\gamma\rho^2(x,y)/s}\frac{\gamma\rho^2(x,y)}{s}\,\mu(dx) \notag \\
&\quad-\int_M |\nabla p_s|^2 \,\e^{\gamma\rho^2(x,y)/s}\rho(x,y)\Delta \rho(\cdot, y)(x)\,\mu(dx).
\end{align}

By the Laplacian comparison theorem (see e.g. <cit.>),
we have

\begin{align*}
-\Delta \rho(\cdot,y)(x)\leq \frac{d-1}{\rho(x,y)}+K\rho(x,y)
\end{align*}

outside of the cut-locus. This yields in particular

\begin{align*}
&-\frac{\gamma}{s}\int_M |\nabla p_s|^2 \,\e^{\gamma\rho^2(x,y)/s}\rho(x,y)\Delta \rho(\cdot,y)(x)\,\mu(dx)\\
&\leq \frac{ (d-1)\gamma}{s}\int_M |\nabla p_s|^2 \,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)+K\int_M |\nabla p_s|^2 \,\e^{\gamma\rho^2(x,y)/s}\,\frac{\gamma\rho^2(x,y)}{s}\,\mu(dx)\\
&\leq  \frac{ (d-1)\gamma}{s}\int_M |\nabla p_s|^2 \,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)+K\int_M |\nabla p_s|^2 \,\e^{\frac{\gamma' \rho^2(x,y)}{s}}\,\mu(dx)\\
&\le  \frac{C_\gz(1+Ks)}{s^2}E(s,y,\gz')
\end{align*}

with $\gamma<\gamma'<2\alpha$ as in (<ref>). Combining this
estimate with (<ref>), it follows

\begin{align}\label{eqn-esI1}
\mathrm{I}_1\leq C_\gamma\frac{(1+Ks)}{s^2}E(s,y,\gz')
\leq \frac{c(1+Ks)}{s^2\sqrt{V(y,\!\sqrt{s})}}\,\e^{2C' s}.
\end{align}

Altogether (<ref>) through (<ref>), we
conclude that

\begin{align*}
\int_M |\Hess p_s|_\HS^2\,\e^{\gamma\rho^2(x,y)/s}\,\mu(dx)
=\mathrm{I}_1+\mathrm{I}_2+\mathrm{I}_3 \leq c\frac{(1+Ks)}{s^2\sqrt{V(y,\!\sqrt{s})}}\,\e^{2C's}.
\end{align*}

which completes the proof of Proposition <ref> by using the

\begin{equation*}
|\Hess p_s|^2\leq |\Hess p_s|^2_\HS.\qedhere
\end{equation*}

Using Proposition <ref>, we obtain the following
$L^2$-integral estimates for Hessian of heat kernel.

Assume (<ref>) holds.  Fix $\az\in (0,\frac{1}{4})$ as in
Lemma .  There exist $0<\beta<\alpha$ and
$C''>C'>0$ with $C'$ as in Lemma <ref> such that

\begin{align*}
\int_{\rho(x,y)\geq t^{1/2}}|\Hess_xp_s(x,y)|\,\mu(dx)\leq
C\big(1+\sqrt{s}\big)\,\e^{C''s}\,\e^{-\beta t/s}s^{-1}
\end{align*}

for all $y\in M$ and $s,t>0$.

Let $0<\beta<\alpha$. By Cauchy's inequality we obtain

\begin{align*}
&\int_{\rho(x,y)\geq t^{1/2}}|\Hess_xp_s(x,y)|\,\mu(dx)\\
& \leq \lf(\int_M|\Hess_xp_s(x,y)|^2\,\e^{2\beta\,\rho^2(x,y)/s}\,\mu(dx)\r)^{1/2}\lf(\int_{\rho(x,y)\geq t^{1/2}}\e^{-2\beta\,\rho^2(x,y)/s}\,\mu(dx)\r)^{1/2}\\
&\leq \frac{C\,\e^{C's}(1+\sqrt{s})}{s\sqrt{V(y,\!\sqrt{s})}}\sqrt{V(y,\!\sqrt{s})}\,\e^{-\beta t/s}\,\e^{\tilde{C}s}\\
&=\frac{C(1+\sqrt{s})}{s}\,\e^{C's-\beta t/s}\,\e^{\tilde{C}s}\leq \frac{C(1+\sqrt{s})}{s}\,\e^{C''s-\beta t/s}
\end{align*}

where the second inequality follows from Proposition <ref>
and inequality (<ref>). This finishes the proof.

§.§ Stochastic Hessian formulas and pointwise estimates for
Hessian of heat kernel

In this subsection, we establish some pointwise and $L^p$-integral
estimates for the Hessian of the heat kernel when $p>2$. To this end, let
us first introduce some necessary notations.  For fixed $x\in M$, let
$B_t$ be the stochastic anti-development of $X.(x)$ which is a
Brownian motion in $T_xM$.  Let $/\!/_t\colon T_x M\to T_{X_t(x)} M$
be the parallel transport and $Q_t\colon T_{x}M\rightarrow T_{X_t}M$ be the damped parallel transport
defined as the solution to the following pathwise ordinary covariant differential
equation along the trajectories of $X_t$,

\begin{align}\label{eqn-Qt}
DQ_t=-\Ric^{\sharp}Q_t\,dt,\quad Q_0={\rm id}_{T_{x}M}
\end{align}

with $DQ_t=/\!/_t\,d\,/\!/_t^{-1}Q_t$.
For each $w\in T_{x}M$
define an operator-valued process
$W_t(\cdot,w): T_{x}M\rightarrow T_{X_t}M$ by

\begin{align*}
W_t(\cdot, w)=Q_t\int_0^tQ_r^{-1}R(/\!/_r\,dB_r,Q_r(\cdot))Q_r(w)-Q_t\int_0^tQ_r^{-1}(\nabla \Ric^{\sharp}+d^*R)(Q_r(\cdot), Q_r(w))\,dr.
\end{align*}

This means that the process $W_t(\cdot,w)$ is the solution to the
following covariant Itô equation

\begin{equation*}\left\{\begin{aligned}
DW_t(\cdot,w)&=R(/\!/_tdB_t, Q_t(\cdot))Q_t(w)-(d^*R+\nabla
\Ric^{\sharp})(Q_t(\cdot), Q_t(w))\,dt -\Ric^{\sharp}(W_t(\cdot,
w))\,dt,\\ W_0(\cdot,w)&=0.
\end{aligned}\right.\end{equation*}

First, we collect some easy curvature estimates for Riemannian
manifolds $M$ satisfying condition (<ref>).

Assume that (<ref>) holds. There exist constants $C>0$,
$\theta>0$ such that for any $t>0$,

\begin{align}\label{R-esti}
\sup_{x\in M}\E^x\left[\int_0^t|R|^2(X_s)\,ds+\int_0^t|\nabla \Ric^{\sharp}+d^*R|^2(X_s)\,ds\right]\leq C\e^{2 \theta t}.
\end{align}

By means of the trivial inequality $s\leq\e^s$ for $s\geq0$, we see

\begin{align*}
&\sup_{x\in M}\E^x\left[\int_0^t|R|^2(X_s)\,ds+\int_0^t|\nabla \Ric^{\sharp}+d^*R|^2(X_s)\,ds\right]\\
&\leq \sup_{x\in M}\E^x\left[\int_0^tK_1(X_s)\,ds+\int_0^tK_2(X_s)\,ds\right]\\
&\leq \sup_{x\in M}\E^x\left[\exp\left(\int_0^tK_1(X_s)\,ds \right)\right]
+\sup_{x\in M}\E^x\left[\exp\left(\int_0^tK_2(X_s)\,ds\right)\right].
\end{align*}

Recall in <cit.> that for any
$K\in \mathcal{K}(M)$, there exist positive constants $c$ and $C$

\begin{align*}
\sup_{x\in M}\E^x\left[\exp\left(\int_0^t |K(X_s)|\,ds\right)\right]\leq C\e^{ct}.
\end{align*}

This immediately completes the proof by substituting $K_1$ and $K_2$
as $K$.

By this lemma and Bismut-type Hessian formula, we obtain the following
Hessian estimates for heat semigroup.

Assume that (<ref>) holds.  Then, for any $p\geq 2$, there
exist positive constants $C$ and $\theta$ as in Lemma
such that, for every $f\in \mathcal{B}_b(M)$,

\begin{align}
&t\,|\Hess P_tf|(x)\leq C(1+\sqrt{t})\,\e^{(2K+\theta)t}\lf(P_t|f|^p\r)^{1/p}(x); \label{lem1-i1} \\
&t\,\big\| |\Hess P_tf|\,\big\|_p\leq C(1+\sqrt{t})\,\e^{(2K+\theta)t}\|f\|_p .\label{lem1-i2}
\end{align}

We start by recalling the following Bismut type formula for
$\Hess P_t$ (see [2, 19, 9]). Let $x\in D$ with
$v,w\in T_xM$, $f\in \mathcal{B}_b(M)$ and $0<s<t$.  Suppose that
$D_1$ and $D_2$ are regular domains such that $x\in D_1$ and
$\bar{D}_1\subset D_2\subset D$, and denote by $\sigma,\tau$ the
exit times of $X.(x)$ from $D_1$ and $D_2$ respectively.  Assume
that $k,\ell$ are bounded adapted processes with paths in the
Cameron-Martin space $L^2([0,t]; [0,1])$ such that $k_r=0$ for
$r\geq \sigma \wedge s$, $k_0=1$, $\ell_r=1$ for
$r\leq \sigma \wedge s$, $\ell_s=0$ for $r\geq \tau \wedge t$.
Then, for $f\in \mathcal{B}_b(M)$, we have

\begin{align}
(\Hess P_tf))(v,w)
&=-\E^x\left[f(X_t)\int_0^t\langle W_s(\dot{k}_sv,  w), /\!/_s \,d B_s\rangle\right]\notag\\
&\quad+\E\left[f(X_t)\int_s^t\langle Q_r(\dot{\ell}_rw), /\!/_r\,  d B_r\rangle
\int_0^{s}\langle Q_r(\dot{k}_rv), /\!/_r\,d B_r\rangle\right].\label{HessBismutFormula}
\end{align}

Note that under condition $\bf (H)$ one has $|Q_t|^2\leq \e^{2Kt}$
and according to the Itô formula, then for $t<\tau<\tau_D$ and

\begin{align*}
d |W_s(v, w)|^2&=2\,\langle R(/\!/_s d B_s, Q_s(v))Q_s(w), W_s(v,w)\rangle+2 |R^{\sharp,\sharp}(Q_s(v),Q_s(w))|_{\rm HS}^2\, d s\\
&\quad- 2\,\langle (d^* R+\nabla \Ric^{\sharp})(Q_s(v), Q_s(w)), W_s(v,w) \rangle\, d s\\
&\quad- 2\Ric(W_s(v,w),W_s(v,w))\,d s\\
&\leq 2\,\langle R(/\!/_s d B_s, Q_s(v))Q_s(w), W_s(v,w)\rangle+ 2\e^{4Ks}K_1(X_s)\, d s\\
&\quad+ 2\sqrt{K_2(X_s)} \,\e^{2Ks} |W_s(v,w)|\, d s+2K |W_s(v,w)|^2\, d s\\
&\leq 2\,\langle R(/\!/_s d B_s, Q_s(v))Q_s(w), W_s(v,w)\rangle+ 2\e^{4Ks}K_1(X_s)\, d s\\
&\quad+ (2\delta)^{-1}\e^{4Ks}K_2(X_s) \, d s+(2K+\delta) |W_s(v,w)|^2\, d s.
\end{align*}

From this and Lemma <ref> there exist constants $C>0$ and
$\theta>0$ such that

\begin{align}\label{est-Wt}
\E|W_{t\wedge \tau_D}(v,w)|^2
&\leq 2\e^{(2K+\delta)t}\left(\E\left[\int_0^t\e^{(2K-\delta)s}K_1(X_s)\,ds\right]+(4\delta)^{-1}\E\left[\int_0^t\e^{(2K-\delta)s}K_2(X_s)\,ds\right]\right) \\
&\leq  2\max\{1,(4\delta)^{-1}\}\e^{4Kt} \E\left[\int_0^tK_1(X_s)\,ds+\int_0^tK_2(X_s)\,ds\right] \notag\\
&\leq  C \,\e^{(4K+2\theta)t}.\notag
\end{align}

Since the upper of $ \E|W_{t\wedge \tau_D}(v,w)|^2$ is independent
of the domain $D$, taking a sequence of domains $D_n$ increasing to
$M$, we conclude that
$ \E|W_{t}(v,w)|^2\leq C \,\e^{(4K+2\theta)t}.$ All these estimates
allow to verify that formula (<ref>) holds for the
deterministic functions $k_s=\frac{t-2s}{t}\vee 0$ and
$\ell_s=1\wedge \frac{2(t-s)}{t}$ as well.  Then, by estimate
(<ref>) and $|Q_s|\leq \e^{Ks}$, we get for any $p\geq 2$,

\begin{align}\label{Hess-pointwise-esti}
|\Hess P_tf|
&\leq  \frac{2}{t}(P_t|f|^p)^{1/p}\left[\int_0^{t/2}\E|W_s(v,w)|^2\,ds\right]^{1/2}\notag\\
&\quad+ \frac{4}{t^2}(P_t|f|^p)^{1/p}\E\left[\Big(\int_0^{t/2}\langle Q_s(w),/\!/_sdB_s\rangle\Big)^2
\E\Big[\Big(\int_{t/2}^t\langle Q_s(v), /\!/_sdB_s \rangle\Big)^2\Big|\SF_{t/2}\Big]\right]^{1/2}\notag\\
&\leq \frac{2}{t}(P_t|f|^p)^{1/p}\left(\int_0^{t/2}\e^{(4K+2\theta)s}\,ds\right)^{1/2}\notag\\
&\quad+ \frac{4}{t^2}(P_t|f|^p)^{1/p}\E\def\ovz{\overrightarrow}
\left[\int_0^{t/2}\e^{2Ks}\,ds\right]^{1/2}\left(\int_{t/2}^t\e^{2Ks}\,ds\right)^{1/2}\notag\\
&\leq  \frac{C}{\sqrt{t}}\,\e^{(2K+\theta)t}(P_t|f|^p)^{1/p}+
\frac{C}{t}\,\e^{2Kt}(P_t|f|^p)^{1/p}\notag\\
&\leq  \frac{C(1+\sqrt{t})}{t}\,\e^{(2K+\theta)t}(P_t|f|^p)^{1/p}
\end{align}

which proves (<ref>).  Moreover, integration with respect to
the measure $\mu$ yields

\begin{align*}
\mu\big(|\Hess P_tf|^p\big)^{1/p}\leq  \frac{C(1+\sqrt{t})}{t}\,\e^{(2K+\theta)t} \mu(|f|^p)^{1/p}.
\end{align*}

This proves (<ref>) and finishes the proof.

The following proposition gives a pointwise estimate for the Hessian
of heat semigroup.

Assume that (<ref>) holds.  Then there exist constants
$C, \theta>0$ such that

\begin{align}\label{hess-gra}
|\Hess P_tf|\leq
\e^{2Kt}\left(P_t|\Hess f|^2\right)^{1/2}+C\e^{(2K+\theta) t}\,\left(P_t|df|^2\right)^{1/2}
\end{align}

for $t>0$ and $f\in C_0^{\infty}(M)$.

Recall the following relation

\begin{align}\label{eqn-Re1}
d\Delta f=\tr \Hess(df)-df(\Ric^{\sharp})
\end{align}

where $\Hess(df)=\nabla^2df$.
Moreover, using the differential Bianchi identity
<cit.>, we have

\begin{align}\label{eqn-Re2}
\Hess(\Delta f)=\tr \Hess(\Hess f)-2(\Hess f)
(\Ric^{\sharp}\otimes {\rm id}-R^{\#,\#})-df(d^*R+\nabla \Ric^{\sharp}).
\end{align}

Now, for $s>0$ and $v,w\in T_xM$, let

\begin{align*}
N_s(v,w):=(\Hess P_{t-s}f)(Q_s(v), Q_s(w))(X_s)+(dP_{t-s}f)(W_s(v,w))(X_s).
\end{align*}

We claim that $N_s$ is a local martingale.  Combining
(<ref>), (<ref>) and applying Itô's formula, we

\begin{align*}
&=\big(\nabla_{/\!/_sdB_s}\Hess P_{t-s}f\big)(Q_s(v), Q_s(w))+(\Hess P_{t-s}f)(\frac{D}{ds}Q_s(v),Q_s(w))\,ds\\
&\quad+(\Hess P_{t-s}f)(Q_s(v),\frac{D}{ds}Q_s(w))\,ds+\partial_s(\Hess P_{t-s}f)(Q_s(v),Q_s(w))\,ds
\\
&\quad+\tr\Hess(\Hess P_{t-s}f)(Q_s(v), Q_s(w))\,ds+(\nabla_{/\!/_sdB_s}d P_{t-s}f)(W_s(v,w))\\
&\quad+(dP_{t-s}f)(DW_s(v,w))+\langle d(dP_{t-s}f), DW_s(v,w) \rangle+\partial_s(dP_{t-s}f)(W_s(v,w))\,ds\\
&\quad+\tr\Hess (dP_{t-s}f)(W_s(v,w))\,ds\mequal 0,
\end{align*}

where $\mequal$ denotes equality modulo the differential of a local
martingale, and here we used for the quadratic covariation the

\begin{align*} [d(dP_{t-s}f), DW_s(v,w)]_s=(\Hess
\end{align*}

Therefore, $N_s$ is a local martingale. Recall that condition
$\bf (H)$ implies

\begin{align*}
|Q_t|^2\leq \e^{2Kt} \quad\text{and}\quad \E|W_t(v,w)|^2\leq C \e^{(4K+2\theta)t}  \end{align*}

for some positive constants. In view of estimate
(<ref>), $|\Hess P_{t-s}f|$ is bounded on
$[0,t-\varepsilon]\times M$, and

$$|\nabla P_{t-s}f|\leq \e^{K(t-s)}\,\|\,|\nabla f|\,\|_{\infty}$$

as a consequence of the derivative formula (see e.g. <cit.>):
for all $v\in T_xM$ and $x\in M$,

\begin{align}\label{derivative-formula}
\langle (\nabla P_tf)(x),v\rangle=\E^x\left[\langle(\nabla f)(X_t),Q_tv\rangle\right].
\end{align}

Thus $N_s$ is a martingale on the time
interval $[0,t-\varepsilon]$, and by taking expectation at time $0$
and $t-\varepsilon$, we arrive at

\begin{align}\label{Hessian-formula}
(\Hess P_tf)(v,w)=\E\Big[(\Hess P_{\varepsilon} f)(Q_{t-\varepsilon}(v),Q_{t-\varepsilon}(w))(X_t)+ dP_{\varepsilon}f(W_{t-\varepsilon}(v,w))(X_{t-\varepsilon})\Big].
\end{align}

Letting $\varepsilon$ tend to $0$, it then follows that

\begin{align*}
|\Hess P_tf(v,w)|
& \leq \Big|\E\Big[\Hess f(Q_t(v), Q_t(w))(X_t)\Big]\Big|+ \Big|\E\Big[df(W_t(v,w))(X_t)\Big]\Big|\\
&\le P_t\Big[|\Hess f|\,|Q_t(v)|\,|Q_t(w)|\Big]+(P_t|df|^2)^{1/2}\E\Big[|W_t(v,w)|^2\Big]^{1/2}\\
&\le (P_t|\Hess f|^2)^{1/2}\,\e^{2Kt}+C\,\e^{(2K+\theta)t}\,(P_t|df|^2)^{1/2}.\qedhere
\end{align*}

$$\|R\|_{\infty}:=\sup_{x\in M} |R|(x)\quad \mbox{ and} \quad
\|\nabla \Ric^{\sharp}+d^*R\|_{\infty}:=\sup_{x\in M} |\nabla
\Ric^{\sharp}+d^*R|(x).$$
If assumption (<ref>) is replaced by
(<ref>) with $K\geq 0$, $\|R\|_{\infty}<\infty$ and
$\|\nabla \Ric^{\sharp}+d^*R\|_{\infty}<\infty$, then it is
straightforward to see from (<ref>) that for $0<\delta<2K$,

\begin{align*}
\E|W_t(v,w)|^2 \leq  \e^{(2K+\delta)t} \E\left[\int_0^t\e^{(2K-\delta)s}\|R\|_{\infty}^2\,ds+\int_0^t\e^{(2K-\delta)s}\|\nabla \Ric^{\sharp}+d^*R\|_{\infty}^2\,ds\right]\leq C\e^{4Kt} t
\end{align*}

for some explicit constant $C>0$. Proceeding as in the proof of
Proposition <ref> above, we obtain in this case

\begin{align*}
|\Hess P_tf|\leq  \e^{2Kt}\big(P_t|\Hess f|^2\big)^{1/2}+C\sqrt{t}\e^{2K t}\big(P_t|df|^2\big)^{1/2}.
\end{align*}

Proposition <ref> is used to establish pointwise Gaussian upper
bounds for the Hessian of the heat kernel via inequality
(<ref>). Note that Carron [7] applied pointwise
Gaussian upper bounds for the gradient of the heat kernel to the Riesz
transform on complete manifolds whose Ricci curvature satisfies a
quadratic decay control. Here, we shall use the following Hessian heat
kernel estimate in the next section to establish the
Calderón-Zygmund inequality.

Assume that (<ref>) holds. Fix $\az\in (0,1/4)$ as in Lemma
and let $\theta>0$ be as in
(<ref>).  Then there exist $0<\beta<2\alpha$ and $C_3>0$ such

\begin{align*}
|\Hess_xp_{t}(x,y)|&\leq \frac{C(1+\sqrt{t})}{tV(y,\!\sqrt{t})}\exp\left(-\beta\,\frac{\rho^2(x,y)}{t}+\frac{1}{2}\big(C_{3}+\theta\big)t\right).
\end{align*}

for any $t>0$ and $x$, $y\in M$.

For $t>0$ and $x$, $y\in M$, write

\begin{align*}
\Hess_xp_{2t}(x,y)=\Hess P_t(p_t(\cdot,y))(x)
\end{align*}

Applying Proposition <ref> and taking $\gz\in (0,2\az)$, we

\begin{align*}
\big|\Hess P_t(p_t(\cdot,y))(x)\big|
&\leq  \e^{2Kt}\left(\int |\Hess_zp_t(z,y)|^2\,\e^{\gamma{\rho^2(z,y)}/t}
\,\e^{-\gamma{\rho^2(z,y)}/t}\,|p_t(x,z)|\,\mu(dz)\right)^{1/2} \\
&\quad+C\e^{(2K+\theta)t}\left(\int |{\nabla}_z p_t(z,y) |^2\,\e^{\gamma{\rho^2(z,y)}/t}\,\e^{-\gamma{\rho^2(z,y)}/t}\,|p_t(x,z)|\,\mu(dz)\right)^{1/2}. \notag
\end{align*}

Moreover, by Lemma <ref> and Proposition <ref>, we
know that there exists $C'>0$ such that

\begin{align*}
\int_M \lf(|{\nabla}_zp_t(z,y)|^2\,\e^{\gamma{\rho^2(z,y)}/t}+ |\Hess_zp_t(z,y)|^2
\,\e^{\gamma{\rho^2(z,y)}/t}\r)\,\mu(dz)\leq \frac{C(1+Kt)}{t^2V(y,\!\sqrt{t})}\,\e^{2C't}.
\end{align*}

This further implies that there exists a constant $C_{3}>0$ such

\begin{align*}%\label{Hess-pointwise-estimate}
\big|\Hess P_t(p_t(\cdot,y))(x)\big|
\frac{C(1+\sqrt{t})}{t\sqrt{V(y,\!\sqrt{t})}}\,\e^{((C'+2K)+\theta)t}\sup_{z\in M}
\left\{\e^{-\gamma{\rho^2(z,y)}/{t}}|p_t(x,z)|\right\}^{1/2}\\
&\le \frac{C(1+\sqrt{t})}{{t}{V(y,\!\sqrt{t})}}\,\e^{(C_{3}+\theta)t}
\,\e^{-\beta {\rho^2(x,y)}/t}
\end{align*}

for $\beta$ small enough, where in the last inequality we have used
the estimate

\begin{align*}
\sup_{z\in M}\Big\{\e^{-\gamma{\rho^2(z,y)}/t}|p_t(x,z)|\Big\}\leq \frac{\e^{C_1Kt}}{V(y,\!\sqrt{t})}\e^{-2\beta {\rho^2(x,y)}/t},
\end{align*}

which can be deduced from Lemma <ref>. Thus, we conclude

\begin{align*}
|\Hess_xp_{2t}(x,y)|\leq \frac{C(1+\sqrt{t})}{tV(y,\!\sqrt{t})}\,\e^{(C_{3}+\theta)t}\,\e^{-\beta{\rho^2(x,y)}/t}
\end{align*}

which finishes the proof of Proposition <ref>.

By interpolating the $L^p$-Hessian inequality with the $L^2$-Gaffney
off-diagonal estimates, one obtains $L^p$-Gaffney off-diagonal
estimates for any $p>2$ as follows.

Assume that (<ref>) hold. Then for $p>2$ there exists
constants $C,C_4>0$ such that for all $t\in (0,\fz)$, all Borel
subsets $E,F\subset M$ with compact closure, and all $f\in L^p(M)$
with supp$f\subset E$,

\begin{align*}
\big\|\mathbbm{1}_Ft |\Hess P_tf|\big\|_p\leq
\e^{-{C_4\rho^2(E,F)}/t}\,\|f\|_p.
\end{align*}

Let $p>2$ and $t>0$. For $E,F$ and $f$ as above, by inequality
(<ref>), there exist positive constants $C$ and $C_2$ such

\begin{align*}
\int_Ft^p|\Hess P_tf|^p(x)\,\mu(dx)
&\leq C\e^{(2K+\theta)pt}(1+\sqrt{t})^p\int_F(P_t|f|^{p/2})^2(x)\,\mu(dx)\\
&\leq C\e^{-{2C_2 \rho^2(E,F)}/t}\,\e^{(2K+\theta)pt}\,(1+\sqrt{t})^p \int_{E} |f|^p\,\mu(dx)
\end{align*}

where the last inequality follows from the $L^2$-Gaffney
off-diagonal estimates for $P_tf$ (Lemma <ref>),
i.e. the existence of positive constants $C$ and $C_2$ such that

\begin{align*}
\big\|\mathbbm{1}_F |P_tf^{p/2}|\,\big\|_{2}\leq C \,\e^{-{C_2 \rho^2(E,F)}/t}\left\|\mathbbm{1}_E f^{p/2}\right\|_{2}
\end{align*}

for all $f\in L^p(M)$ with $\text{supp} f\subset E$.

§ CALDERÓN-ZYGMUND INEQUALITY FOR $P\NEQ 2$

In this section, we prove the main results of this paper. We first
prove Theorem <ref> in Section <ref>, namely CZ($p$)
for $p\in (1,2)$ under the condition (<ref>); then in Section
<ref>, we prove Theorem <ref>, namely, CZ($p$)
for $p>2$ under the condition (<ref>). Recall that CZ($2$)
always holds under condition (<ref>).

§.§ The case $p\in (1,2)$

Let $M$ be a complete Riemannian manifold satisfying (<ref>).
In this subsection, we prove Theorem <ref> and show that the
Calderón-Zygmund inequality holds for $p\in (1,2)$. To be precise,

\begin{align}\label{eqn-defT}
T=\Hess\,(\Delta+\sz)^{-1}=\dint_0^\fz \e^{-\sz t} \Hess P_t\,dt.
\end{align}

We show that $T$ is bounded on $L^p(M)$ for any $p\in (1,2)$, provided
$\sigma>0$ is large enough. Since $T$ is already bounded on $L^2(M)$,
using interpolation, it suffices to prove that $T$ is of weak type
$(1,1)$, that is for some constant $C$,

\begin{align}\label{eqn-bdT}
\mu\big(\{x\in M\colon |Tf(x)|>\lambda\}\big)\leq C\dfrac{\|f\|_1}{\lambda}
\end{align}

for all $\lambda>0$ and $f\in L^1(M)$. To this end, we need the
following technical lemma from <cit.> on the
finite overlap property of $M$.

Assume that (<ref>) holds. There exists a countable subset
$\mathcal{C}=\{x_j\}_{j\in \Lambda}\subset M$ such that

(i) $M=\cup_{j\in \Lambda} B(x_j,1)$;

(ii) $\{B(x_j,1/2)\}_{j\in\Lambda}$ are disjoint;

(iii) there exists $N_0\in\N$ such that for any
$x\in M$, there are at most $N_0$ balls $B(x_j,4)$ containing $x$;

(iv) for any $c_0\geq 1$, there exists $C>0$ such that
for any $j\in\Lambda$, $x\in B(x_j,c_0)$ and $r\in (0,\fz)$,

\begin{align*}
\mu\lf(B(x,2r)\cap B(x_j,c_0)\r)\le C \mu\lf(B(x,r)\cap B(x_j,c_0)\r)
\end{align*}

\begin{align*}
\mu(B(x,r))\le C \mu\lf(B(x,r)\cap B(x_j,c_0)\r)
\end{align*}

for any $x\in B(x_j,c_0)$ and $r\in (0,2c_0]$.

The following lemma provides the localization argument to prove

Assume that (<ref>) holds.  Let
$\mathcal{C}=\{x_j\}_{j\in\Lambda}$ be a countable subset of $M$
having finite overlap property as in Lemma
. Let $\sigma>C''$ where $C''$ is as in
Corollary <ref>. Suppose that there exists a positive constant
$C$ such that

\begin{align}\label{eqn-localT}
\mu\lf(\{x\colon  \mathbbm{1}_{B(x_j,2)}|Tf(x)|>\lambda\}\r)\leq \frac{C}{\lambda}\|f\|_1
\end{align}

for any $j\in \Lambda$, $\lz\in (0,\fz)$ and
$f\in C_0^\fz (B(x_j,1))$. Then (<ref>) holds for any
$f\in C_0^\fz(M)$.

For $j\in \Lambda$, let $B_j:=B(x_j,1)$ and let
$\{\varphi_j\}_{j\in\Lambda}$ be a $C^{\infty}$-partition of unity
such that $0\le \varphi_j\le 1$ and is supported in $B_j$.  Then,
for any $f\in C_0^{\infty}(M)$ and $x\in M$, write

\begin{align*}
Tf(x)= \sum_{j\in\Lambda} \mathbbm{1}_{2B_j}T(f\varphi_j)(x)+\sum_{j\in\Lambda}(1-\mathbbm{1}_{2B_j})T(f\varphi_{j})(x),
\end{align*}

which yields that for any $\lambda>0$,

\begin{align*}
\mu(\{x\colon |Tf(x)|>\lambda\})
&\leq  \mu\bigg(\bigg\{x\colon \dsum_{j\in\Lambda} \mathbbm{1}_{2B_j}|T(f\varphi_j)(x)|>\frac{\lambda}{2}\bigg\}\bigg)
+ \mu\bigg(\bigg\{x\colon \dsum_{j\in\Lambda}(1-\mathbbm{1}_{2B_j})|T(f\varphi_j)(x)|>\frac{\lambda}{2}\bigg\}\bigg)\\
\end{align*}

For $\mathrm{I}_1$, by Lemma <ref>(iii) and
(<ref>), we have

\begin{align}\label{lem-eI1}
\mathrm{I}_1\le  \dsum_{j\in\Lambda}\mu\bigg(\bigg\{x:  \mathbbm{1}_{2B_j}|T(f\varphi_j)(x)|>\frac{\lambda}{2N_0}\bigg\}\bigg)
\ls\frac{1}{\lambda}\|f\|_1
\end{align}

as desired, where the notation $a\ls b$ means $a\leq Cb$ for some
constant $C$.

To bound $\mathrm{I}_2$, again by Lemma <ref>(iii), since
$\varphi_j$ is supported in $B_j$, it is easy to see that

\begin{align*}
\dsum_{j\in\Lambda}|(1-\mathbbm{1}_{2B_j})(x)\varphi_j(y)|\leq N_0 \mathbbm{1}_{\{\rho(x,y)\geq 1\}}.
\end{align*}

Hence, according to the definition of $T$ in (<ref>) and
Corollary <ref>, we get

\begin{align*}
\mathrm{I}_2
&\leq \frac{2}{\lz} \dsum_{j\in\Lambda} \lf\|\,\big|\lf(1-\mathbbm{1}_{2B_j}\r)
\\
&\ls \frac{1}{\lz} \int_M\lf(\int_0^{\infty}\e^{-\sigma t} \int_M \lf|\Hess_xp_t(x,y)\r|\sum_{j\in\Lambda}|(1-\mathbbm{1}_{2B_j})(x)\varphi_j(y)||f(y)|\,\mu(dy)\,dt\r)\mu(dx)\\
&\ls\frac{1}{\lz}\int_{M}\int_0^{\infty}\e^{-\sigma t} \lf(\int_{\rho(x,y)\geq 1} |\Hess_xp_t(x,y)|\,\mu(dx)\r)\,dt\, |f(y)|\,\mu(dy)\\
&\leq \frac{1}{\lz} \int_M |f(y)|\,\mu(dy) \int_0^{\infty}\e^{-\sigma t}(1+\sqrt{t})\,\e^{C''t}\,\e^{-\beta/ t}t^{-1}\,dt,
\end{align*}

where $\bz\in (0,\az)$. Thus, since $\sigma >C''$, we obtain

\begin{align*}
\mathrm{I}_2\ls \frac{1}{\lz}\int_0^{\infty}\e^{(C''-\sigma )t
-{\beta}/{t}}\frac{(1+\sqrt{t})}{t}\,dt \,\|f\|_1\ls \frac{1}{\lz}\|f\|_1,
\end{align*}

which combined with the estimate about $I_1$ in (<ref>)
finishes the proof of Lemma <ref>.

To prove property (<ref>), we remove the subscript $j$ and
simply write $B$ for each $B_j:=B(x_j,1)$.  Let $c_0\geq 1$. By Lemma
<ref>(iv), we have that $(c_0B,\mu,\rho)$ is a metric measure
subspace satisfying the volume doubling property that there
exists $C_D\ge 1$ such that

\begin{align}\label{eqn-D}
\mu\lf(B(x,2r)\cap c_0B\r)\leq C_D \,\mu\lf(B(x,r)\cap c_0B\r) \tag{\bf D}
\end{align}

for all $x\in c_0B$ and $r>0$.

We also need the following Calderón-Zygmund decomposition from

Let $(\mathcal{X},\nu,\rho)$ be a metric measure space satisfying
(<ref>) with $c_0B$ replaced by $\mathcal{X}$.  Let
$f\in L^1(\mathcal{X})$ and $\lz\in (0,\fz)$. Assume
$\|f\|_{L^1}< \lz \nu(\mathcal{X})$.  Then $f$ has a decomposition
of the form
$$\textstyle f=g+b=g+\sum_ib_i$$
such that

* $g(x)\leq C\lambda$ for almost all $x\in M$;

* there exists a sequence of balls $B_i=B(x_i,r_i)$ so that the
support of each $b_i$ is contained in $B_i$:

\begin{align*}
\int_{\mathcal{X}} |b_i(x)|\, \nu(dx)\leq C\lambda \nu(B_i)\quad \text{and}\quad
\int_{\mathcal{X}} b_i(x)\, \nu(dx)=0;
\end{align*}

* $\displaystyle\sum_i \nu(B_i)\leq
\frac{C}{\lambda}\int_{\mathcal{X}} |f(x)|\, \nu(dx)$;

* there exists $k_0\in \mathbb{N}^*$ such that each point of $M$
is contained in at most $k_0$ balls $B_i$.

Assume that (<ref>) holds.  Let $\lz\in (0,\fz)$ and
$f\in L^1(B)$ be as in Lemma . Assume that
$\{b_i\}$ is the sequence of bad functions as in Lemma
and $\{P_{t}^{\sz}\}_{t\ge 0}$ the heat
semigroup associated to $-(\Delta+\sz)$ with $\sz>0$. Then there
exist $\sigma >0$ and $C>0$ independent of $f$ such that

\begin{align*}
\Big\|\sum_iP^{\sigma}_{t_i}b_i\Big\|_2^2\leq C\lambda \|f\|_1
\end{align*}

where $t_i=r_i^2$ with $r_i$ denoting the radius of the ball $B_i$
as in Lemma <ref>(b).

Recall that $\supp b_i \subset B(x_i, \sqrt{t_i})$. Using the upper
bound of the heat kernel in Lemma <ref> and
Lemma <ref>, we have for $x\in M$,

\begin{align*}
&\leq \int_M \frac{\e^{-\sz' t_i-\alpha \frac{\rho^2(x,y)}{t_i}}}{V(x,\sqrt{t_i})}|b_i(y)|\,\mu(dy)\\
& \leq  \frac{C}{V(x,\!\sqrt{t_i})} \,\e^{-\sz' t_i-\alpha'\frac{\rho^2(x,x_i)}{t_i}}\int_{B_i} |b_i(y)|\,
\mu(dy)\\
&\leq  C_2 \lambda \int_M \frac{\e^{-\sz' t_i-\alpha{''}\frac{{ \rho^2(x,y)}}{t_i}}}{V(x,\sqrt{t_i})}\mathbbm{1}_{B_i}(y)\,\mu(dy),
\end{align*}

where $\sz'=\sz-C_1K>0$ and $0<\az''<\az'<\az$.  It is therefore
sufficient to verify that

\begin{align}\label{lem-le1}
\biggl\|\sum_i \int_M \frac{\e^{-\sz' t_i-\alpha'' \frac{\rho^2(\cdot, y)}{t_i}}}
{V(\cdot, \sqrt{t_i})}\mathbbm{1}_{B_i}(y)\,\mu(dy)\biggr\|_2\ls\biggl\|\sum_i\mathbbm{1}_{B_i}\biggr\|_2,
\end{align}

since from this and Lemma <ref> we obtain as
consequence, as desired,

\begin{align*}
\Big\|\sum_iP_{t_i}^{\sigma}b_i\Big\|^2_2\ls\lambda^2\biggl\|\sum_i\mathbbm{1}_{B_i}\biggr\|^2_2\ls \lambda^2 \sum_i \mu(B_i)\ls \lambda \|f\|_1.
\end{align*}

In order to prove (<ref>), we write by duality

\begin{align}\label{lem-le2}
&\biggl\|\sum_i \int_M \frac{\e^{-\sz' t_i-\alpha''\frac{\rho^2(\cdot,y)}{t_i}}}{V(\cdot,\!\sqrt{t_i})}\mathbbm{1}_{B_i}(y)\,\mu(dy)\biggr\|_2\\ \notag
&\qquad =\sup_{\|u\|_2=1}\lf|\int_M\lf(\sum_i\int_M\frac{\e^{-\sz' t_i-\alpha'' \frac{\rho^2(x,y)}{t_i}}}{V(x,\sqrt{t_i})}\mathbbm{1}_{B_i}(y)\,\mu(dy)\r)u(x)\,\mu(dx)\r|\\ \notag
\leq \sup_{\|u\|_2=1}\int_M \sum_i\lf(\int_M \frac{\e^{-\sz' t_i-\alpha''\frac{\rho^2(x,y)}{t_i}}}{V(x,\sqrt{t_i})}|u(x)|\,\mu(dx)\r)\mathbbm{1}_{B_i}(y)\,\mu(dy).
\end{align}

By (<ref>), we have for any $x\in M$ and $y\in B_i$,

\begin{align*}
V(y,\!\sqrt{t_i})\leq C\lf(1+\frac{\rho(x,y)}{\sqrt{t_i}}\r)^d \e^{C(t_i^{1/2}+\rho(x,y))}
\end{align*}

with $\dz\in [0,2)$. From this, we obtain that there exist
$0<\tilde{\az}<\az'''<\az''$ such that

\begin{align*}
&\int_M \frac{\e^{-\sz' t_i-\alpha'' \frac{\rho^2(x,y)}{t_i}}}{V(x,\!\sqrt{t_i})}\,|u(x)|\,\mu(dx)\\
&\hs\ls\frac{\e^{-\frac{1}{2}\sz' t_i}}{V(y,\!\sqrt{t_i})}\int_M \e^{-\alpha'''\frac{\rho^2(x,y)}{t_i}}|u(x)|\,\mu(dx)\\
+\sum_{k=0}^\fz\int_{2^k\!\sqrt{t_i}\leq \rho(x,y)< 2^{k+1}\!\sqrt{t_i}}\e^{-\alpha'''\frac{\rho^2(x,y)}{t_i}}|u(x)|\,\mu(dx)\Bigg)\\
&\hs\leq \frac{1}{V(y,\!\sqrt{t_i})}\lf(\int_{B(y,\!\sqrt{t_i})}|u(x)|\,\mu(dx)+\sum_{k=0}^\fz\e^{-\alpha''' 2^{2k}}\int_{B(y,2^{k+1}\sqrt{t_i})}|u(x)|\,\mu(dx)\r)\\
&\hs=\lf(1+\sum_{k=0}^\fz\frac{V(y,\!2^{k+1}\sqrt{t_i})}{V(y,\!\sqrt{t_i})}\,\e^{-\alpha'''2^{2k}}\r)(\SM u)(y)\\
&\hs\leq \left(1+C\sum_{k=0}^\fz2^{(k+1)d}\e^{C2^{k+1}\sqrt{t_i}}\e^{-\tilde{\az} 2^{2k}}\right)(\SM u)(y)\\
& \hs\leq \left(1+C\sum_{k=0}^\fz2^{(k+1)d}\e^{C2^{k+1}}\e^{-\tilde{\az} 2^{2k}}\right)(\SM u)(y)\ls (\SM u)(y),
\end{align*}

\begin{align*} (\SM
\end{align*}

denotes the Hardy-Littlewood maximal function of $u$.  This
together with (<ref>) and the $L^2$-boundedness of
$\SM $ gives

\begin{align*}
\biggl\|\sum_i\int_M \frac{\e^{-\alpha''\frac{d^2(\cdot, y)}{t_i}}}{V(\cdot,\!\sqrt{t_i})}\mathbbm{1}_{B_i}(y)\,\mu(dy)\biggr\|_2\ls
\sup_{\|u\|_2=1}\int_M (\SM u)(y)\sum_{i}\mathbbm{1}_{B_i}(y)\,\mu(dy)
\ls\biggl\|\sum_i\mathbbm{1}_{B_i}\biggr\|_2,
\end{align*}

which shows that (<ref>) holds true and finishes the
proof of Lemma <ref>.

With the help of Lemmas <ref> through <ref>, we
are now in position to the proof of Theorem <ref>.

Recall that $T=\Hess(\Delta+\sz)^{-1}$. We choose $\sigma$ big
enough such that $\sigma>\max\{C'', C_1K\}$, where $\tilde{\alpha}$
is as in proof of Lemma <ref>.  By Lemma <ref>, it
suffices to prove

\begin{align}\label{eqn-localT-1}
\mu\big(\{x\in 2B\colon  |Tf(x)|>\lambda\}\big)\ls \dfrac{\|f\|_1 }{\lambda}
\end{align}

for all $f\in C_0^\fz(B)$. By means of Lemma <ref> with
$\mathcal{X}=B$, we deduce that $f$ has a decomposition

$$\textstyle f=g+b=g+\sum_ib_i$$
which implies

\begin{align}\label{eqn-them-8}
\mu\big(\{x\in 2B: |Tf(x)|>\lambda\}\big)
&\leq \mu\lf(\lf\{x\in 2B:|Tg(x)|>\frac{\lambda}{2}\r\}\r)+
\mu\lf(\lf\{x\in 2B: |Tb(x)|>\frac{\lambda}{2}\r\}\r)\\ \notag
\end{align}

Using the facts that $T$ is bounded on $L^2(M)$ and that
$|g(x)|\leq C\lambda$, we obtain as desired

\begin{align}\label{eqn-them-9}
\mathrm{I}_1\ls\lambda^{-2}\|Tg\|_2^2\ls\lambda^{-2}\|g\|_2^2
\ls\lambda^{-1}\|g\|_1\ls\lambda^{-1}\|f\|_1.
\end{align}

We now turn to the estimate of $\mathrm{I}_2$. Recall that
$\{P_{t}^{\sz}\}_{t\geq 0}$ is the heat semigroup generated by
$-(\Delta+\sz)$, that is $P_t^{\sz}=\e^{-t\sz}P_t$.  We write

\begin{align*}
\end{align*}

where $t_i=r_i^2$ with $r_i$ the radius of $B_i$.  By Lemma
<ref>, we have

\begin{align*}
\biggl\|\sum_iP_{t_i}^{\sz}b_i\biggr\|_2^2\ls\lambda \|f\|_1.
\end{align*}

This combined with the $L^2$-boundedness of $T$ yields

\begin{align}\label{eqn-them1-1}
\mu\left(\left\{x\in 2B\colon \lf|T\lf(\sum_i P_{t_i}^{\sz}b_i\r)(x)\r|>\frac{\lambda}{2}\right\}\right)
\ls \frac{1}{\lambda}\|f\|_1
\end{align}

as desired. Consider now the term $T\sum_i(I-P_{t_i}^{\sz})b_i$. We

\begin{align}\label{eqn-them-4}
&\mu\lf(\lf\{x\in 2B\colon\lf|T\lf(\sum_i(I-P_{t_i}^{\sz})b_i\r)(x)\r|>\frac{\lambda}2\r\}\r)\\ \notag
&\hs\leq \sum_{i}\mu(2B_i)+\mu\lf(\lf\{x\in 2B\setminus\cup_i2B_i\colon \lf|T\lf(\sum_i(I-P_{t_i}^{\sz})b_i\r)\r|(x)>\frac{\lambda}{2}\r\}\r).
\end{align}

From Lemma <ref>, it follows that

\begin{align}\label{eqn-them1-2}
\sum_{i}\mu(2B_i)\ls\frac{\|f\|_1}{\lambda}
\end{align}

as desired. To estimate the second term, let $k_{t_i}^{\sz}(x,y)$
denote the integral kernel of the operator $T(I-P_{t_i}^{\sz})$.
Note that

\begin{align*}
\end{align*}

\begin{align*}
T(I-P_{t_i}^{\sz})&=\Hess(\Delta+\sigma)^{-1}(I-P_{t_i}^{\sz})=\int_0^{t_i}\Hess P_{s}^{\sz}\,ds.
\end{align*}

\begin{align}\label{eqn-kernelKt}
\end{align}

where $p_s^{\sigma}=\e^{-\sz s}p_s$ is the kernel of $P_s^{\sz}$
with respect to $\mu$.  Since $b_i$ is supported in $B_i$, we have

\begin{align}\label{eqn-them-3}
\int_{2B\setminus (2B_i)}
\lf|T\lf((I-P_{t_i}^{\sz})b_i\r)(x)\r|\,\mu(dx)
\int_{2B\setminus (2B_i)}\lf(\int_{B_i}|k_{t_i}^{\sz}(x,y)||b_i(y)|\,\mu(dy)\r)\,\mu(dx)\\ \notag
&\leq \int_{B_i}\lf(\int_{\rho(x,y)\geq t_i^{1/2}}|k_{t_i}^{\sz}(x,y)|\,\mu(dx)\r)|b_i(y)|\,\mu(dy).
\end{align}

Now by means of (<ref>) and Corollary <ref>, we

\begin{align*}
\int_{\rho(x,y)\geq t_i^{1/2}}|k_{t_i}^{\sz}(x,y)|\,\mu(dx)
\int_0^{t_i}\lf(\int_{\rho(x,y)\geq t_i^{1/2}}|\Hess_xp_{s}(x,y)|\,\mu(dx)\r)\e^{-s\sigma}\,ds\\
&\leq C\int_0^{t_i} \e^{-\beta t_i/s}s^{-1}\,\e^{C''s}\,\e^{-s\sigma}(1+\sqrt{s})\,ds\\
&\leq C\int_0^1\frac{\e^{-\beta/u}}{u}\,du<\fz,
\end{align*}

where for the last inequality we use the fact that

$$\e^{s(C''-\sigma) }(1+\sqrt{s})<\fz,\quad s\in (0,\fz).$$

The estimate above together with (<ref>) and Lemma
<ref> implies that

\begin{align}\label{eqn=them-6}
\mu\lf(\lf\{x\in 2B\setminus\cup_i2B_i\colon
\lf|T\lf(\sum_i(I-P_{t_i}^{\sz})b_i\r)(x)\r|>\frac{\lambda}{2}\r\}\r)
\ls \frac{\|f\|_1}{\lz}.
\end{align}

Altogether, combining (<ref>) through (<ref>),
(<ref>) and (<ref>), we conclude that
(<ref>) holds which completes the proof of Theorem

§.§ The case $p\in (2,\infty)$

Let $M$ be a complete Riemannian manifold satisfying (<ref>).
In this subsection, we prove Theorem <ref> and show that the
Calderón-Zygmund inequality CZ$(p)$ holds for all
$p\in(2,\fz)$, that is,

\begin{align}\label{eqn-Tpg2}
\lf\|T f\r\|_{p}\ls \|f\|_p
\end{align}

holds for any $f\in L^p(M)$ with $T$ as in (<ref>) and
$\sigma\in (2K+\theta,\infty)$ where $K,\theta$ are respectively in
(<ref>) and (<ref>).

To this end, let $w$ be a $C^{\infty}$ function on $[0,\infty)$
satisfying $0\leq w\leq 1$ and

\begin{align*}
w(t)=\begin{cases}1\ & \text{on} \  [0,3/4],\\
0 \ & \text{on} \ [1,\fz),
\end{cases}
\end{align*}

and let $\wz{T}$ be an operator defined by

\begin{align}\label{eqn-wzT}
\wz{T}f:=\int_0^{\infty}v(t)\Hess P_tf\,dt
\end{align}

with $v(t):=w(t)\e^{-\sigma t}$. We need the following lemma, which
reduces (<ref>) to a time and spatial localized version.

Assume that (<ref>) holds.  Let $p\in (2,\infty)$ and
$\{x_j\}_{j\in \Lambda}$ be a countable subset of $M$ having the
finite overlap property as in Lemma . If there
exists a positive constant $C$ such that

\begin{align}\label{eqn-wTfb}
\lf\| \, |\wz T (f)|\, \r\|_{L^p(B(x_j,4))}\le C \lf\|f\r\|_{L^p(B(x_j,1))}
\end{align}

holds for any $j\in \Lambda$ and $f\in C_0^\fz(B(x_j,1))$ with
$\wz T$ defined as in (<ref>), then (<ref>) holds.

By the fact that $w\equiv1$ on $[0,3/4]$, we obtain from Proposition
<ref> that if $\sigma> 2K+\theta$, then for any $g\in L^p(M)$,

\begin{align*}
\Big\|\int_0^{\infty}(1-w(t))\,\e^{-\sigma t}|\Hess P_{t}g|\,dt\Big\|_{p}\ls
\int_{3/4}^{\infty}\e^{(2K+\theta-\sigma)t}\, \frac{1+\sqrt{t}}{t}dt\, \|g\|_{p}\ls \|g\|_p.
\end{align*}

This and (<ref>) imply that to prove (<ref>), it
suffices to show that

\begin{align}\label{eqn-wTfb-1}
\lf\| \,|\wz T (g)|\, \r\|_{p}\ls \lf\|g\r\|_{p}
\end{align}

for any $g\in C^\fz_0(M)$.

Let $(x_j)_{j\in \Lambda}$ be a countable subset of $M$ having the
finite overlap property as in Lemma <ref>.  Let $\{\fai\}_j$
be a corresponding $C^{\infty}$ partition of unity such that
$0\leq \varphi_j\leq 1$ and $\fai_j$ is supported in
$B_j:=B(x_j,1)$. Let $\chi_j$ be the characteristic function of the
ball $4B_j$. For any $g\in C_0^{\infty}(M)$ and $x\in M$, write

\begin{align}\label{eqn-wztq}
\wz{T}g(x)\leq \sum_{j\in \Lambda} \chi_j\wz{T}(g\varphi_j)(x)+\sum_{j\in \Lambda}(1-\chi_j)\wz{T}(g\varphi_j)(x)=:{\rm I}(x)+{\rm II}(x).
\end{align}

By Lemma <ref>, we know

\begin{align*}
\sum_{j\in \Lambda}\lf|(1-\chi_j)(x)\varphi_j(y)\r|\leq N_0 \mathbbm{1}_{\{\rho(x,y)\geq 3\}}.
\end{align*}

Hence, by Hölder's inequality, we have

\begin{align*}
{\rm II}(x)
&\leq \int_0^1 \int_M \lf|\Hess_xp_t(x,y)\r|\lf(\sum_{j\in \Lambda}\lf|(1-\chi_j)(x)\varphi_j(y)\r|\r)|g(y)|\,\mu(dy)\,dt\\
&\leq N_0\int_0^1 \int_{\rho(x,y)\geq 3}\lf|\Hess_xp_t(x,y)||g(y)\r|\,\mu(dy)\,dt\\
\left(V(y,\!\sqrt{t})\right)^{p/p'}|g(y)|^p\,\mu(dy)\right)^{1/p}\frac{\e^{-{c}/{t}}}{t}\,
\end{align*}

where $c={\gamma p'}/{p}$ and $\gz$ is a positive constant small
enough so that in view of Proposition <ref> and
an argument similar to the proof of (<ref>), it implies
that for any $t\in (0,1)$,

\begin{align*}
\int_M |t\Hess_xp_t(x,y)|^p\,\e^{\gamma\frac{\rho^2(x,y)}{t}}\,\mu(dx)\ls \frac{1}{\left(V(y,\!\sqrt{t})\right)^{p-1}}.
\end{align*}

This immediately implies

\begin{align}\label{eqn-wztqp}
&\int_M |{\rm II}(x)|^p\,\mu(dx)\\ \notag
&\ls\int_M \left[\int_0^1\left(\int_M |t\Hess_x p_t(x,y)|^p\,\e^{\gamma{\rho^2(x,y)}/{t}}V(y,\!\sqrt{t})^{p/p'}|f(y)|^p\,\mu(dy)\right)^{1/p}
\frac{\e^{-c/t}}{t}\,dt\right]^p\,\mu(dx)\\ \notag
&\ls\int_M \int_0^1 \left(\int_M |t\Hess_x p_t(x,y)|^p\,\e^{\gamma{\rho^2(x,y)}/{t}}\,V(y,\!\sqrt{t})^{p/p'}|f(y)|^p\,\mu(dy)\right)\,dt\,\mu(dx)\\ \notag
&\qquad\times\left(\int_0^1\frac{\e^{-{cp'}/{t}}}{t^{p'}}\,dt\right)^{p/p'}\\ \notag
&\ls \int_0^1\lf[\int_M[V(y,\!\sqrt{t})]^{p/p'}|f(y)|^p\lf(\int_M |t\Hess_x p_t(x,y)|^p\,\e^{\gamma{\rho^2(x,y)}/{t}}\,\mu(dx)\r)\,\mu(dy)\r]\,dt \\ \notag
&\ls \int_M |f(y)|^p\mu(dy)
\end{align}

as desired.

Now we estimate $\mathrm{I}(x)$. Using Lemma <ref>, we know
that the balls $\{4B_j\}_{j\in \Lambda}$ are of uniform overlap and

\begin{align*}
\sum_j \|h\chi_j\|_{p'}^{p'}\ls \|h\|_{p'}^{p'}
\end{align*}

for all $h\in C_0^\fz(M)$. Since $g\fai_j\in C_0^\fz (B(x_j,1))$ and
using (<ref>), we conclude that

\begin{align*}
\lf|\int_Mh(x)\mathrm{I}(x)\,\mu(dx)\r|
\int_M|h(x)|\Big|\sum_j\chi_j\wz{T}(g\varphi_j)(x)\Big|\,\mu(dx)\\
&\ls \sum_j \|g\varphi_j\|_p\|h\chi_j\|_{p'}\ls \|g\|_{p}\|h\|_{p'},
\end{align*}

which together with (<ref>) and (<ref>) implies
(<ref>), and hence finishes the proof of Lemma

From Lemma <ref>, it is easy to see that to prove $T$ is of
strong type $(p,p)$ for $p>2$, it suffices to show that $\wz T$ is
bounded from $L^p(B_j,\mu)$ to $L^p(4B_j,\mu)$ for each $j\in \Lambda$
as in (<ref>). To this end, we need an $L^p$ local bounded
criterion from [3] via maximal functions.  Recall that
the local maximal function  by

\begin{align}\label{eqn-localMF}
(\SM_{\rm loc}f)(x):=\sup_{\substack{B\ni x \\
r(B)\leq 32}}\frac{1}{\mu(B)}\int_B |f|\,d\mu,\quad x\in M,
\end{align}

for any locally integrable function $f$ on $M$.  From (<ref>),
it follows that $\SM_{\rm loc}$ is bounded on $L^p(M)$ for all
$1<p\leq \infty$.  For a measurable subset $E\subset M$, the maximal function relative to $E$  is defined by

\begin{align}\label{eqn-relaMF}
(\SM_Ef)(x):=\sup_{B\text{ ball in }  M,\,B\ni x} \frac{1}{\mu(B\cap E)}\int_{B\cap E}|f|\,d\mu,
\quad x\in E,
\end{align}

for any locally integrable function $f$ on $M$. If in particular $E$
is a ball with radius $r$, it is enough to consider balls $B$ with
radii not exceeding $2r$. It is also easy to see $\SM_E$ is of weak
type $(1,1)$ and $L^p(M)$-bounded for $1<p\leq \infty$ if $E$
satisfies the relative doubling property, namely, if there
exists a constant $C_E$ such that for all $x\in E$ and $r>0$,

\begin{align}\label{RD}
\mu(B(x,2r)\cap E)\leq C_E\,\mu(B(x,r)\cap E).
\end{align}

Note that in Lemma <ref> (iv), for any $j\in \Lambda$, the
subsets $4B_j$ satisfy the relative doubling property (<ref>) with
a constant independent of $j$.

The following theorem is essential to the proof of Theorem

Let $p\in (2,\infty)$ and assume that (<ref>) holds. For any
ball $B$ of $M$ centered in $4B_j$ with radius less than 8, assume

* there exists an integer $n$ depending only on condition
(<ref>) such that the map
$f\rightarrow \SM^{\#}_{4B_j,\tilde{T},n}f$ is bounded from
$L^p(B_j,\mu)$ to $L^p(4B_j,\mu)$ with operator norm independent
of $j$, where

\begin{align*}
\SM^{\#}_{4B_j,\tilde{T},n}f(x):=\dsup_{B \text{ \rm ball in }M,\, B\ni x}\lf(\frac{1}{\mu(B\cap 4B_j)}\dint_{B\cap 4B_j}
\lf|\wz{T}(I-P_{r^2})^nf(y)\r|^2\,\mu(dy)\r)^{1/2};
\end{align*}

* for all $k\in \{1,2,\ldots,n\}$, and all $f\in L^2(M,\mu)$
supported in $B_j$, there exists a sublinear operator $S_j$
bounded from $L^p(B_j,\mu)$ to $L^p(4B_j,\mu)$ with operator norm
independent of $j$ such that for $x\in B\cap 4B_j$,

\begin{align}\label{eq2}
\left(\frac{1}{\mu(B\cap 4B_j)}\int_{B\cap 4B_j}|\wz{T}P_{kr^2}f|^p\,d\mu\right)^{1/p}
\leq C\left(\SM_{4B_j}(|\wz{T}f|^2)+(S_jf)^2\right)^{1/2}(x)
\end{align}

where $\SM_{4B_j}$ is as in (<ref>).

Then $\wz T$ is bounded from $L^p(B_j,\mu)$ to $L^p(4B_j,\mu)$, that
is (<ref>) holds with a constant depending on $p$, the
doubling constant $C_D$ in (<ref>), the operator norm of
$\wz{T}$ on $L^2(M,\mu)$, the operator norms of
$\SM^{\#}_{4B_j,\tilde{T},n}$ and $S_j$ on $L^p$, and the constant
in (<ref>).

We use the following result <cit.>:

$p_0\in (2,\infty]$ and $(M,\mu,\rho)$ be a measured metric
space. Suppose that $T$ is a bounded sublinear operator which is
bounded on $L^2(M,\mu)$, and let $\{A_r\}_{r>0}$, be a family of
linear operators acting on $L^2(M,\mu)$. Let $E_1$ and $E_2$ be two
subsets of $M$ such that $E_2$ has the relative doubling property,
$\mu(E_2)<\infty$ and $E_1\subset E_2$. Assume that

* the sharp maximal functional $ \SM_{E_2,T,A}^{\#}$
is bounded from $L^p(E_1,\mu)$ into $L^p(E_2,\mu)$ for all
$p\in (2,p_0)$, where

\begin{align*}
&(\SM^{\#}_{E_2,T,A}f)^2(x)=\sup_{B \ \text{ball in}\ M,\ B\ni x}\frac{1}{\mu(B\cap E_2)}\int_{B\cap E_2}|T(I-A_{r(B)})f|^2\,d\mu
\end{align*}

for $x\in E_2$;

* for some sublinear operator $S$ bounded from $L^p(E_1,\mu)$
into $L^p(E_2, \mu)$ for all $p\in (2,p_0)$,

\begin{align}\label{TA-2}
\left(\frac{1}{\mu(B\cap E_2)}\int_{B\cap E_2}|TA_{r(B)}f|^{p_0}\,d\mu\right)^{1/{p_0}}
\leq C\left(\SM_{E_2}(|Tf|^2)+(S(f)^{2}\right)^{1/2}(x),
\end{align}

for all $f\in L^2(M,\mu)$ supported in $E_1$, all balls $B$ in $M$
and all $x\in B\cap E_2$, where $r(B)$ is the radius of $B$.

If $2<p<p_0$ and $Tf\in L^p(E_2,\mu)$ whenever $f\in L^p(E_1,\mu)$,
then $T$ is bounded from $L^p(E_1,\mu)$ into $L^p(E_2,\mu)$ and its
operator norm is bounded by a constant depending only on the
operator norm of $T$ on $L^2(M,\mu)$, $C_{E_2}$ (see in (<ref>)
for $E_2$), $p$ and $p_0$, the operator norms of
$\SM^{\#}_{4B_j,\tilde{T},n}$ and $S$ on $L^p$, and the constant in

Here we may take $E_1$ and $E_2$ as $B_j$ and $4B_j$ respectively, as
the sets $B_j$ and $4B_j$ have the relatively volume doubling property
as in (<ref>) with the constant $C_E$ independent of $j$ (see Lemma
<ref>). Moreover, taking the operators $\{A_r\}_{r>0}$ such

\begin{align*}
\end{align*}

for some integer $n$ sufficiently large, then the result follows

Now it suffices to check (i) and (ii) of Lemma <ref>. We
establish two technical lemmas which verify (i) and (ii)
respectively. To this end, observe that (<ref>) implies: for
all $r_0>0$ there exists $C_{r_0}$ such that for all $x\in M$,
$r\in (0,r_0)$,

\begin{align*}
V(x,2r)\leq C_{r_0} V(x,r).
\end{align*}

An easy consequence of the definition is that for all $y\in M$,
$0<r<8$ and $s\geq 1$ satisfying $sr<32$,

\begin{align}\label{local-doubling}
V(y,sr)\leq Cs^{D_{L}} V(y,r),
\end{align}

for some constants $C$ and $D_{L}>0$. The following lemma plays an
important role, when checking (i) of Lemma <ref>.

Assume that (<ref>) hold.  Then there exists an integer $n$
such that the inequality

\begin{align}\label{eq-Hess-1}
\dsup_{B \, ball\ in\  M,\,B\ni x}\lf(\frac{1}{\mu(B\cap 4B_j)}\dint_{B\cap 4B_j}
\lf|\wz{T}(I-P_{r^2})^nf(y)\r|^2\,\mu(dy)\r)^{1/2}
\leq C \left(\SM_{\rm loc}(|f|^2)(x)\right)^{1/2}
\end{align}

holds for any $x\in 4B_j$, $f\in L^2(4B_j)$ satisfying
$\supp f\subset 4B_j$, where $\SM_{\rm loc}$ is defined by

Viewing the left-hand side of (<ref>) as the maximal
function relative to $4B_j$, since the radius of $4B_j$ is $4$, it
is enough to consider balls $B$ of radii not exceeding $8$.  Let
$B=B(x_0,r)$ be an arbitrary ball containing $x$ and satisfying
$x_0\in 4B_j$ and $r\in (0,8)$. By Lemma <ref>, we know that

\begin{align}\label{B-control}
\mu(B)\ls \mu(B\cap 4B_j)
\end{align}

with an implicit constant independent of $B$ and $j$. Hence, it is
easy to see

\begin{align*}
\frac{1}{\mu(B\cap 4B_j)}\int_{B\cap 4B_j} \big|\wz{T}(I-P_{r^2})^nf\big|^2\,d\mu\ls \frac{1}{\mu(B)}\int_{B}\big|\wz{T}(I-P_{r^2})^nf\big|^2\,d\mu.
\end{align*}

Thus, we only need to show that

\begin{align}\label{eqn-lem-term1-1}
\dsup_{B\ \text{ball in}\ M,\, B\ni x}\lf(\frac{1}{\mu(B)}\dint_{B}
\lf|\wz{T}(I-P_{r^2})^nf(y)\r|^2\,\mu(dy)\r)^{1/2}
\ls  \left(\SM_{\rm loc}(|f|^2)(x)\right)^{1/2},
\end{align}

for any $x\in 4B_j$.  Since $r<8$, we choose $i_r\in \Z_+$

\begin{align}\label{def-ir}
2^{i_r}r\leq 8< 2^{i_r+1}r.
\end{align}

Denote by

\begin{align}\label{C-B}
C_i\ \text{ the annulus}\  2^{i+1}B\setminus (2^iB)\ \  \ \text{if}\  i\geq 2\, \text{
and} \, C_1=4B.
\end{align}

Using the fact $\supp f\subset 4B_j\subset 2^i B$ when $i>i_r$, we
find that

\begin{align*}
f=\dsum_{i=1}^{i_r} f\mathbbm{1}_{C_i}=:\dsum_{i=1}^{i_r} f_i
\end{align*}

which then implies

\begin{align}\label{eqn-lem-term1-2}
\lf\|\,| \wz{T}(I-P_{r^2})^nf|\,\r\|_{L^2(B)}\leq \dsum_{i=1}^{i_r}\lf\|\,|\wz{T}(I-P_{r^2})^nf_i|\,\r\|_{L^2(B)}.
\end{align}

For $i=1$ we use the $L^2$-boundedness of $\wz{T}(I-P_{r^2})^n$ to

\begin{align}\label{eqn-lem-term1-3}
\lf\|\,|\wz{T}(I-P_{r^2})^nf_1|\,\r\|_{L^2(B)}\leq \|f\|_{L^2(4B)}\leq \mu(4B)^{1/2}(\SM_{\rm loc}(|f|^2)(x))^{1/2}
\end{align}

as desired.  For $i\geq 2$, we infer from (<ref>) that

\begin{align*}
\wz{T}(I-P_{r^2})^nf_i
&=\dint_0^\fz v(t)\dsum_{k=0}^n \binom nk (-1)^k \Hess\,P_{t+kr^2}(f_i)\,dt \\
&=\dint_0^\fz \lf(\dsum_{k=0}^n \binom nk (-1)^k \mathbf{1}_{\{t>kr^2\}} v(t-kr^2)\r)\Hess\, P_{t}(f_i)\,dt\\
&=:\int_0^{\infty}g_r(t)\,\Hess P_tf_i\,dt.
\end{align*}

For $g_r$, according to the definition $v(t)=w(t)\e^{-\sigma t}$ and
an elementary calculation, we observe that

\begin{align*}
\begin{cases}
|g_r(t)|\ls 1 \quad  &\mbox{for}\ \  0<t \leq \ (1+nr^2)\wedge (1+n)r^2,\\
|g_r(t)| \ls r^{2n} \quad  &\mbox{for}\ \   (1+nr^2)\wedge (1+n)r^2<t \leq 1+nr^2,\\
g_r(t)=0\quad & \mbox{for}\ \ t > 1+nr^2,
\end{cases}
\end{align*}

where the second estimate comes from $w^{(i)}(t)\ls 1$ for
$t\in [0,1]$ and $i\in \{1,\ldots,n\}$, along with

\begin{align*}
\Bigg|\sum_{k=0}^n\binom nk (-1)^k w(t-kr^2) \Bigg|\leq C_n \sup_{u\geq \frac{t}{n+1}}\Big|(w(u)\,\e^{-\sigma u})^{(n)}\Big|r^{2n}\ls r^{2n}
\end{align*}

for all $t\geq (1+n)r^2$ and some constant $C_n$.  Combined with
Lemma <ref>, this gives

\begin{align*}
\lf\|\,\big|\wz{T}(I-P_{r^2})^nf_i\big|\,\r\|_{L^2(B)}\ls \left(\int_0^{\infty}|g_r(t)|(1+\sqrt{t})\,\e^{-{\alpha'4^ir^2}/{t}}\,\frac{ dt}{t}\right)\|f_i\|_{L^2(C_i)}
\end{align*}

where using the fact that $0<r<8$, we know

\begin{align*}
&\ls  \int_{0}^{(1+nr^2)\wedge (1+n)r^2}(1+\sqrt{t})\,\e^{-\frac{\alpha' 4^i r^2}{t}}\frac{dt}{t}+\int_{(1+nr^2)\wedge (1+n)r^2}^{1+nr^2}(1+\sqrt{t})r^{2n}\,\e^{-\frac{\alpha' 4^i r^2}{t}}\frac{dt}{t} \\
&\leq C_n \left(\int_0^{(n+1)r^2}\e^{-\frac{\alpha' 4^i r^2}{t}}\frac{dt}{t} +\int_{(1+nr^2)\wedge (1+n)r^2}^{1+nr^2}(1+\sqrt{t})r^{2n}\frac{t^{n-1}}{4^{in}r^{2n}}\, {dt}\right)\\
&\leq C_n\left(4^{-in}+ 4^{-in}r^{2n}(1+\sqrt{r})\right)\leq C_n'4^{-in}.
\end{align*}

Now an easy consequence of the local doubling
(<ref>), since $r(2^iB)\leq 8$ when
$1\leq i\leq i_r$, is that

\begin{align*}
\mu(2^{i+1}B)\leq C2^{(i+1)D_{L}}\mu(B),
\end{align*}

with constants $C$ and $D_{L}$ independent of $B$ and
$i$. Therefore, as $C_i\subset 2^{i+1}B$,

\begin{align*}
\|f\|_{L^2(C_i)}\leq \mu(2^{i+1}B)^{1/2}(\SM_{\rm loc}(|f|^2)(x))^{1/2}\leq C2^{iD_{L}/2}\mu(B)^{1/2}(\SM_{\rm loc}(|f|^2)(x))^{1/2}.
\end{align*}

Using the definition of $i_r$, $r\leq 8$, and by choosing
$2n>D_{L}/2$, we finally obtain

\begin{align*}
&\big\|\,|\wz{T}(I-P_{r^2})^nf|\,\big\|_{L^2(B)}\leq C'\left(\sum_{i=1}^{i_r}2^{i(D_{L}/2-2n)}\right)\mu(B)^{1/2}(\SM_{\rm loc}(|f|^2)(x))^{1/2},
\end{align*}

which proves Proposition <ref>.

The following lemma is essential to the proof of part (ii) of Lemma

Let $p\in (2,\infty)$. Assume that (<ref>) holds. For a ball
$B$ with radius $r\in (0,8)$, let $i=i_r$ be an integer such that
(<ref>) holds. Then the following estimate hold: For
any $C^2$-function $f$ supported in $C_i$ as in
(<ref>), and each $k\in \{1,\ldots, n\}$, where
$n\in \mathbb{N}$ is chosen according to Lemma
, one has

\begin{align}
&\left(\frac{1}{\mu(B)}\int_B|\Hess P_{kr^2}f|^p\,d\mu\right)^{1/p}\label{Hess-f-esti1}\\
&\quad\leq C\e^{-\alpha_1 4^i}\left[\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\nabla f|^2\,d\mu\right)^{1/2}+\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\Hess f|^2\,d\mu\right)^{1/2}\right]\notag  \\
\intertext{and}
&\left(\frac{1}{\mu(B)}\int_B|\Hess P_{kr^2}f|^p\,d\mu\right)^{1/p}\leq   \frac{C\e^{-\alpha_2 4^i}}{r^2} \left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}| f|^2\,d\mu\right)^{1/2}  \label{Hess-f-esti2}
\end{align}

for some positive constants $C$ and $\alpha_1$, $\alpha_2$ depending
on $K$,$\theta$, $p$ and the constants in (<ref>).

We first observe that Lemma <ref> yields

\begin{align}\label{Hess-ineq}
\left(\int_B|\Hess P_{t}f|^p\,d\mu\right)^{1/p}
& \leq \e^{2Kt}\left(\int_B \Big(P_t|\Hess f|^2\Big)^{p/2}(x)\,\mu(dx)\right)^{1/p}\notag\\
&\quad +C \e^{(2K+\theta)t}\left(\int_B  \left(P_t|df|^2\right)^{p/2}(x)\,\mu(dx)\right)^{1/p}.
\end{align}

We substitute $t=kr^2$ in estimate (<ref>). By Lemma
<ref>, one has the upper bound of $p_t(x,y)$,

\begin{align*}
p_t(x,y)\leq \frac{C_{\alpha}}{V(y,\!\sqrt{t})}\exp{\left(-\alpha\frac{\rho^2(x,y)}{t}+C_1Kt\right)}
\end{align*}

for all $x,y\in M$. Since $r\leq 8$, it follows from the above
estimate that for all $x\in B$,

\begin{align*}
\left|P_{kr^2}(|df|^2)(x)\right|
&\leq C\int_{C_i}V(y,\sqrt{k}r)^{-1}\exp\left(-\alpha\frac{\rho^2(x,y)}{kr^2}+C_1Kkr^2\right)|df|^2(y)\,\mu(dy)\\
&\leq C\e^{-{\alpha}4^i/k}\int_{C_i}V(y,\sqrt{k}r)^{-1}|df|^2(y)\,\mu(dy).
\end{align*}

Moreover, since $y\in C_i$, we have
$2^{i+1}B\subset B(y, 2^{i+2}r)$, and then by (<ref>), we
know that

\begin{align*}
\frac{1}{V(y,\sqrt{k}r)}\leq \frac{2^{d(i+2)}\e^{C2^{i+2}}}{V(y,2^{i+2}r)}\leq \frac{2^{d(i+2)}\e^{C2^{i+2}}}{\mu(2^{i+1}B)}.
\end{align*}

It then follows that

\begin{align*}
|P_{kr^2}|df|^2(x)|\leq C\e^{-c4^i}2^{d(i+2)}\e^{C2^{i+2}}\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\nabla f|^2\,d\mu\right)
\end{align*}

for all $x\in B$, and there exists $\alpha_1<c$ such that

\begin{align}\label{gradient-ineq-esti}
\left(\frac{1}{\mu(B)}\int_{B}(P_{kr^2}|\nabla f|^2)^{p/2}\,d\mu\right)^{1/p}\leq C\e^{-\alpha_1 4^i}\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\nabla f|^2\,d\mu\right)^{1/2}.
\end{align}

With similar arguments, we obtain

\begin{align*}
\left(\frac{1}{\mu(B)}\int_{B}(P_{kr^2}|\Hess f|^2)^{p/2}\,d\mu\right)^{1/p}\leq C\e^{-\alpha_1 4^i}\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\Hess f|^2\,d\mu\right)^{1/2}.
\end{align*}

Altogether these yields

\begin{align*}
&\left(\frac{1}{\mu(B)}\int_B|\Hess P_{kr^2}f|^p\,d\mu\right)^{1/p}\\
&\leq C \e^{-\alpha_1 4^i}\left(\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\nabla f|^2\,d\mu\right)^{1/2}+\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\Hess f|^2\,d\mu\right)^{1/2}\right),
\end{align*}

which completes the proof of (<ref>).

Next we observe from Proposition <ref> that

\begin{align}\label{Hess-ineq2}
\left(\int_B|\Hess P_{t}f|^p\,d\mu\right)^{1/p}
&\leq\frac{(1+\sqrt{t})}{t}\e^{(2K+\theta)t}\left(\int_B \Big(P_t|f|^2\Big)^{p/2}(x)\,\mu(dx)\right)^{1/p}.
\end{align}

Substituting $t=kr^2$ in (<ref>), as $r\in (0,8)$,
there exists a constant $C$ such that

\begin{align*}
\left(\int_B|\Hess P_{kr^2}f|^p\,d\mu\right)^{1/p}
& \leq \frac{C}{r^2}\left(\int_B \Big(P_{kr^2}|f|^2\Big)^{p/2}(x)\,\mu(dx)\right)^{1/p}.
\end{align*}

With similar arguments as for (<ref>), we then
complete the proof of (<ref>).

For any $p\in (2,p_0)$. Fix such a $p$. By expanding
$I-(I-P_{r^2})^n$ it suffices to show

\begin{align*}
\left(\frac{1}{\mu(B)}\int_B |\Hess P_{kr^2}f|^p\right)^{1/p}\leq C\SM (|\Hess f|^2)^{1/2}(y)
\end{align*}

for $f$ with $f$, $\nabla f$ locally square integrable, $B$ any ball
with $r=r(B)$, $y\in B$ and $k=1,2,\ldots, n$.  Recall that $n$ is
chosen larger than $D_{L}/4$ where $D_{L}$ is given.

On the other hand,

\begin{align*}
&P_t|\nabla f|^2(x)\leq C\SM (|\nabla f|^2)(x), \ \mbox{and}\ \
P_t|\Hess f|^2(x)\leq C\SM (|\Hess f|^2)(x)
\end{align*}

by means of combining the upper estimate

\begin{align*}
p_t(x,y)\leq \frac{C}{V(y,\!\sqrt{t})}\exp\Bigg(-c\frac{\rho^2(x,y)}{t}\Bigg)
\end{align*}

and [?,Proposition 2.4]

Let $(x_j)_{j\in J}$ be a maximal 1-separated subset of $M$; the
collection of balls $B_j=B(x_j,1)$, $j\in J$ covers $M$, whereas the
balls $B(x_j,1/2)$ are pairwise disjoint. It follows the doubling
property that there exists $N\in \mathbb{N}^*$ such that every
$x\in M$ is contained in at most $N$ balls $4B_j=B(x_j,4)$.

Consider a $C^{\infty}$ partition of unity $\varphi_j,j\in J$, such
that $\varphi_j\geq 0$ and is supported in $B_j$. Let
$\mathbbm{1}_j$ be the characteristic function of the ball
$4B_j$. For $f\in C_0^{\infty}(M)$ and $x\in M$, write

\begin{align*}
\wz{T}f(x)\leq \sum_j\mathbbm{1}_j\wz{T}(f\varphi_j)(x)+\sum_j(1-\mathbbm{1}_j)\wz{T}(f\varphi_j)(x)={\rm I}+{\rm II}.
\end{align*}

Let us first treat the term ${\rm II}$. The first observation is
that, from the finite overlap property of the balls $B_j$, we have

\begin{align*}
\sum_j|(1-\mathbbm{1}_j)(x)\varphi_j(y)|\leq N\mathbbm{1}_{\rho(x,y)\geq 3}.
\end{align*}

\begin{align*}
{\rm II}&\leq \int_0^1\int_M |\Hessp_t(x,y)|\sum_j |(1-\mathbbm{1}_j)(x)\varphi_j(y)||f(y)|\,\mu(dy)\,dt\\
&\leq N \int_0^1 \int_{\rho(x,y)\geq
\end{align*}

Form there, we follow the argument by inserting the Gaussian terms
and using the Hölder inequality to estimate the integral on
$\rho(x,y)\geq 3$. Since (LY) for small times applies, we have

\begin{align*}
\int_{\rho(x,y)\geq 3}\frac{\e^{-\frac{p'\gamma\,\rho^2(x,y)}{pt}}}{V(y,\!\sqrt{t})}\,\mu(dy)\leq C\e^{-\frac{c}{t}}
\end{align*}

for some constants $C,c>0$ independent of $t$ and $x$. Therefore,

\begin{align*} {\rm II} \leq C\int_0^1\Big(\int_M|t\Hess _x
\frac{\rho^2(x,y)}{t}}[V(y,\!\sqrt{t})]^{p/p'}|f(y)|^p\,\mu(dy)\Big)^{1/p}
\frac{\e^{-\frac{c}{t}}}{t}\,dt.
\end{align*}

It follows from (Hess$_{p_0}$) that for some $\gamma >0$ and all
$t\in (0,1)$ and $y\in M$,

\begin{align*}
\int_M |\Hess_x p_t(x,y)|^p\,\e^{\gamma'\frac{\rho^2(x,y)}{t}}\,\mu(dy)\leq \frac{C}{t^{p}[V(y,\!\sqrt{t})]^{p-1}}.
\end{align*}

Since the measure $\e^{-c/t}/t\,dt$ has finite mass, one can use
Jensen's inequality with respect to $t$, Fubini's theorem and the
weighted $L^p$ estimate above to conclude that
$\int_M |\mathrm{{\rm II}}|^p\,\mu(dx)$ is bounded by $C\|f\|_p$.

We now turn to $I$ which is the main term. The uniform overlap of
the balls $4B_j$ implies

\begin{align*}
\sum_j\|g\mathbbm{1}_j\|_r^r\leq C\|g\|_r^r
\end{align*}

for all $g$ and $1\leq r\leq \infty$. Hence

\begin{align*}
&\int_M |g(x)|\Big|\sum_j \mathbbm{1}_j\wz{T}(f\varphi_j)(x)\Big|\,\mu(dx)\\
&\leq C\sum_j \|f\varphi_j\|_p\|g\|_{p'}\leq C\|f\|_p\|g\|_q,
\end{align*}

provided we show that

\begin{align*}
\|\mathbbm{1}_j|\wz{T}(f\varphi_j)|\|_p\leq C\|f\varphi_j\|_p
\end{align*}

with a bound uniform in $j$. In other words, we want to show that
$\wz{T}$ maps $L^p(B_j)$ into $L^p(4B_j)$.  To this end we apply
Theorem with $E_1=B_j$ and $E_2=4B_j$ since $4B_j$ has the doubling
property by the following lemma.

With the help of the Lemmas <ref>, <ref> and <ref>,
we are now in position to turn to the proof of Theorem <ref>.

By Lemma <ref>, we only need to show (i) and (ii) of
Lemma <ref> hold true under our condition
(<ref>). We first verify (i) of Lemma
<ref>. Observe from Lemma <ref> that there
is an integer $n$ depending only on $D_{L}$ as in
(<ref>) such that for all $f\in L^p(B_j)$ and
$x\in 4B_j$,

\begin{align*}
\sup_{ B \,\text{ball in}\, M, B\ni x} \frac{1}{\mu(B\cap 4B_j)}\int_{B\cap 4B_j} |\wz{T}(I-P_{r^2})^nf|^2\,d\mu\leq C\SM_{\rm loc}(|f|^2)(x).
\end{align*}

Recall that $\SM_{\rm loc}$ is bounded on $L^p(M)$ for
$1<p\leq \infty$; thus for all $p>2$, the operator
$\SM^{\#}_{4B_j,\tilde{T},n}$ is bounded from $L^p(B_j,\mu)$ to
$L^p(4B_j,\mu)$ uniformly in $j$, i.e.  assertion (i) is proved.

Next, we prove (ii) of Lemma <ref>.  Assume that
$f\in L^2(B_j)$ and let $h=\int_0^{\infty}v(t)P_tf\,dt$ with $v$
defined as in (<ref>). Since $\wz{T}f=\Hess\, h$ and the
inequality (<ref>) for $B\cap 4B_j$, we have

\begin{align*}
\left(\frac{1}{\mu(B\cap 4B_j)}\int_{B\cap 4B_j}|\wz T P_{kr^2}f|^p\,d\mu\right)^{1/p}
&= \left(\frac{1}{\mu(B\cap 4B_j)}\int_{B\cap 4B_j}|\Hess P_{kr^2}h|^p\,d\mu\right)^{1/p}\\
&\ls \left(\frac{1}{\mu(B)}\int_B |\Hess P_{kr^2}h|^p\,d\mu\right)^{1/p}.
\end{align*}

We write

\begin{align*}
\Hess P_{kr^2}h=\sum_{i=1}^{\infty}\Hess P_{kr^2}(h\mathbbm{1}_{C_i})=\sum_{i=1}^{\infty}\Hess P_{kr^2}g_i,
\end{align*}

where $g_i=h\mathbbm{1}_{C_i}$ and $C_i$ is as in (<ref>). Next,
we distinguish the two regimes $i\leq i_r$ and $i>i_r$ where $i_r$
is defined as in (<ref>). In the regime $i\leq i_r$, by the
inequality (<ref>) in Lemma <ref>, we

\begin{align*}
&\left(\frac{1}{\mu(B)}\int_B|\Hess P_{kr^2}g_i|^p\,d\mu\right)^{1/p}\\
&\leq C\e^{-\alpha_1 4^i}\left(\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\nabla h|^2\,d\mu\right)^{1/2}+\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\Hess\, h|^2\,d\mu\right)^{1/2}\right)\\
&\leq  C\e^{-\alpha_1 4^i}\left(\left(\SM_{\rm loc}(|\nabla h|^2)(x)\right)^{1/2}+\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|\Hess\, h|^2\,d\mu\right)^{1/2}\right).
\end{align*}

On the other hand, since $C_i\subset 2^{i+1}B$, we know

\begin{align*}
&\frac{1}{\mu(2^{i+1}B)}\int_{C_i}\big|\Hess\, h\big|^2\,d\mu\\
&\leq \frac{1}{\mu(2^{i+1}B\cap 4B_j)}\int_{(2^{i+1}B)\cap 4B_j}|\Hess\, h|^2d\mu+\frac{1}{\mu(2^{i+1}B)}\int_{2^{i+1}B}\mathbbm{1}_{M\setminus 4B_j}|\Hess\, h|^2\,d\mu\\
&\leq \SM_{4B_j}(|\Hess\, h|^2)(x)+\SM_{\rm loc}(|\Hess\, h|^2\mathbbm{1}_{M\setminus 4B_j})(x),
\end{align*}

for any $x\in B\cap 4B_j$. Hence in this case

\begin{align}\label{part-1}
&\left(\frac{1}{\mu(B)}\int_B \big|\sum_{i=1}^{i_r}\Hess P_{kr^2}g_i\big|^p\,d\mu\right)^{1/p} \\
&\leq \sum_{i=1}^{ i_r}C\e^{-\alpha 4^i} \Big(\SM_{\rm loc}(|\nabla h|^2)+\SM_{4B_j}(|\Hess\,h|^2)+\SM_{\rm loc}(|\Hess\,h|^2\mathbbm{1}_{M\setminus 4B_j})\Big)^{1/2}(x).\notag
\end{align}

For the second regime $i>i_r$, we proceed with inequality (<ref>) in Lemma <ref> that there exist
positive constants $c_1$ and $c_2$ such that

\begin{align}\label{Hess-gi}
\left(\frac{1}{\mu(B)}\int_B|\Hess P_{kr^2}g_i|^p\,d\mu\right)^{1/p}\leq \frac{c_1\e^{-c_2 4^i}}{r^2}\left(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|h|^2\,d\mu\right)^{1/2}.
\end{align}

On the other hand, since $i>i_r$, it is easy to see
$4B_j\subset2^{i+1}B$, thus

\begin{align}\label{h-esti}
\Big(\frac{1}{\mu(2^{i+1}B)}\int_{C_i}|h|^2\,d\mu\Big)^{1/2}&\leq \left(\frac{1}{\mu(2^{i+1}B)}\int_0^1 v(t)\int_{C_i} P_t|f|^2\, d\mu \, dt \right)^{1/2}\\
&\leq C\left(\frac{1}{\mu(4B_j)}\int_{B_j}|f|^2\,d\mu\right)^{1/2}=C\left(\frac{1}{\mu(4B_j)}\int_{B_j}|f|^2\,d\mu\right)^{1/2} \notag\\
&\leq C\big(\SM_{4B_j}(|f|^2)(x)\big)^{1/2}. \notag
\end{align}

The contribution of the terms in the second regime $i>i_r$ is
bounded by combining (<ref>) and (<ref>),

\begin{align}\label{part-2}
\sum_{i>i_r}\left(\frac{1}{\mu(B)}\int_B|\Hess P_{kr^2}g_i|^p\,d\mu\right)^{1/p}\leq \sum_{i>i_r} \frac{c_1\e^{-c_2 4^i}}{r^2}\big(\SM_{4B_j}(|f|^2)(x)\big)^{1/2}
\end{align}

and it remains to recall that $1/r^2\leq 4^i/8$ and
$1/r\leq 2^i/{(2\sqrt{2})}$ when $i>i_r$.

We then conclude from (<ref>) and (<ref>) that for any
$p>2$ and $k\in \{1,2,\ldots,n\}$, there exists a constant $C$
independent of $j$ such that

\begin{align*}
\left(\frac{1}{\mu(B\cap 4B_j)}\int_{B\cap 4B_j}|\wz{T}P_{kr^2}f|^p\,d\mu\right)^{1/p}\leq C(\SM_{4B_j}(|\wz{T}f|^2)+(S_jf)^2)^{1/2}(x)
\end{align*}

all $f\in L^2(M,\mu)$ supported in $B_j$, all balls $B$ in $M$ and
all $x\in B\cap 4B_j$, where the radius $r$ of $B$ is less than 8,
and where

\begin{align}\label{eqn-squref}
(S_jf)^2:=\SM_{\rm loc}(|\wz{T}f|^2\mathbbm{1}_{M\setminus 4B_j})
+\SM_{\rm loc}(|\nabla h|^2)+\SM_{4B_j}(|f|^2).
\end{align}

Our last step is to show that the operator $S_j$ defined in
(<ref>) is bounded from $L^p(B_j)$ to $L^p(4B_j)$ for any
$p\in (2,\fz)$ with operator norm independent of $j$.  By
(<ref>), we only need to show that the operators

$$\big(\SM_{\rm loc}(|\wz{T}f|^2\mathbbm{1}_{M\setminus 4B_j})\big)^{1/2},\ \
\big(\SM_{\rm loc}(|\nabla h|^2)\big)^{1/2} \ \mbox{ and }\
\big(\SM_{4B_j}(|f|^2)\big)^{1/2}$$
are respectively bounded from
$L^p(B_j)$ to $L^p(4B_j)$.  Indeed, for any $f\in L^p(4B_j)$, by Lemma
<ref> we know that $4B_j$ satisfies the doubling property
(<ref>), which combined with $p>2$ implies that
$\big(\SM_{4B_j}(|f|^2)\big)^{1/2}$ is bounded from $L^p(B_j)$ to
$L^p(4B_j)$ by a constant depending only on the doubling property
(<ref>).  On the other hand, using the fact from Bismut's
formula [5] that

\begin{align*}
\lf\| \sqrt t\nabla P_t f \r\|_{p} \ls \e^{Kt} \|f\|_p,
\end{align*}

we deduce from $p>2$ and $\sigma>2K+\theta$ that

\begin{align*}
\big\|\lf(\SM_{\rm loc}(|\nabla h|^2)\r)^{1/2}\big\|_p
\ls \|\nabla h\|_{p}\leq \left(\int_0^{\infty}\frac{v(t)\,\e^{ Kt}}{\sqrt{t}}\,dt\right)\,\|f\|_{L^p(B_j)}\leq C\|f\|_{L^p(B_j)}.
\end{align*}

Finally, the $L^p$-boundedness of

$$\big(\SM_{\rm loc}(|\wz{T}f|^2\mathbbm{1}_{M\setminus
follows from the $L^{p/2}$-boundedness of
$\big(\SM_{\rm loc}(|\wz{T}f|^2\big)^{1/2}$ and an argument similar to
(<ref>). This implies that the operator $S_j$ is bounded
from $L^p(B_j)$ to $L^p(4B_j)$ with an upper bound independent of $j$.

We conclude that the requirements (i) and (ii) in Lemma
<ref> both hold under the condition (<ref>).
Thus, the operator $\wz{T}$ is bounded from $L^p(B_j,\mu)$ to
$L^p(4B_j,\mu)$ for $p>2$ with a constant independent of $j$.
Therefore, by Lemma <ref>, the operator $T$ is strong type
$(p,p)$ for $p>2$. This proves Theorem <ref>.

§ OTHER DISCUSSION

For each $p\in (1,\infty)$, we say that $(T_p)$ holds, if the
transform $|\Hess \Delta^{-1} |$ is bounded on $L^p(M)$.

\begin{align*}
\Hess\Delta^{-1}f=\int_0^{\infty}\Hess P_tf\, {d t}.
\end{align*}

One says that the heat kernel satisfies a Gaussian upper bound if
there exists $C,c>0$ such that for all $t>0$ and $x,y\in M$

\begin{equation}
p_{t}(x, y) \leq \frac{C}{V(x,\! \sqrt{t})} \exp \left\{-\frac{d^{2}(x, y)}{c t}\right\} .
\end{equation}

The case $p>2$ is more difficult. Note that if $(T_p)$ holds, then
it follows from the analytic property for the heat semigroup

\begin{align*}
\||\Hess P_t|\|_{p\rightarrow p}\leq \||\Hess\Delta^{-1}\Delta P_t|\|_{p\rightarrow p}\leq \frac{C}{t}
\end{align*}

for all $t>0$.

We now introduce another inequality:there exists $C>0$ such that for
any ball $B$ with radius $r_B$ and any harmonic function $u$ on
$3B$, it holds

\begin{align*}
\left(\fint_B |\Hessu|^p\,d\mu\right)^{1/p}\leq \frac{C}{r_B}\fint_{2B}|u|\,d\mu
\end{align*}

Let $M$ be a complete non-compact Riemannian manifold satisfying (D)
and . If for some $p_0\in (2,\infty]$ there exists $C>0$ such that,
for all $t>0$,

\begin{align*}
\||\Hess P_t|\|_{p_0\rightarrow p_0}\leq \frac{C}{t},
\end{align*}

then $T$ is bounded on $L^p$ for $2<p<p_0$.

Let $M$ be a complete non-compact Riemannian manifold satisfying
(LY). Let $p_0\in (2,\infty]$. The following assertions are

(1) For all $p\in (2,p_0)$, there exists $C_p$ such that

\begin{align*}
\||\Hess_x  P_t|\|_{p\rightarrow p}\leq \frac{C}{t},
\end{align*}

for all $t>0$.

(2) The operator $T$ is bounded on $L^p$ for $p\in (2,p_0)$.

Let $M$ be a complete non-compact Riemannian manifold satisfying (D)
and (DUE).  If there exists $C$ such that, for all $x,y\in M$,

\begin{align*}
|\Hess_x p_t(x,y)|\leq \frac{C}{t[V(y,\!\sqrt{t})]},
\end{align*}

then the operator $T$ is bounded on $L^p$ and the equivalence holds
for $1<p<\infty$.

For all $r_0>0$, for all $x\in M$, $\sigma>1$, $r\leq r_0$,

\begin{align*}
V(x, \sigma r)\leq m(\sigma)V(x,r)
\end{align*}

Let $M$ be a complete non-compact Riemannian manifold satisfying
$(E)$ and $(P_{loc})$.  If for some $p_0\in (2,\infty]$ and
$\alpha\geq 0$, and for all $t>0$,

\begin{align*}
\||\Hess P_t|\|_{p_0\rightarrow p_0}\leq \frac{C\e^{\alpha t}}{t},
\end{align*}

then $\Hess (\Delta+a)^{-1/2}$ is bounded on $L^p$ for $2<p<p_0$ and

Assume that $M$ is a complete non-compact Riemannian manifold
satisfying (FK). Let $2<p_0\leq \infty$. The following assertions
are equivalent:

(1) ($G_p$) holds for all $2<p<p_0$.

(2) For all $2<p<p_0$, for all $y\in M$ and $t>0$,

\begin{align*}
\||\Hess_xp_t(\cdot,y)|\|_p\leq \frac{C_p}{t[V(y,\!\sqrt{t})]^{1-\frac{1}{p}}}.
\end{align*}

Suppose that $\Ric\geq -K^2$. We consider the operators

\begin{align*}
d(-\Delta+K^2)^{-\frac{1}{2}}: L^2(M)\rightarrow L^2\Lambda^1(M),
\end{align*}

\begin{align*}
D(\Delta_1+K^2)^{-\frac{1}{2}}: L^2\Lambda^1(M)\rightarrow L^2T_0^2(M),
\end{align*}

where $\Delta_1=d\delta+\delta d$ is the Hodge Laplacian on 1-forms,
$D$ denotes its covariant derivative and $T_0^2(M)$ is the vector
bundle of 2-covariant tensor fields. Since, obviously,

\begin{align*}
-\Delta+K^2=D^*D+K^2\geq D^*D+K^2
\end{align*}

and, by the Weitzenböck formula,

\begin{align*}
\Delta_1+K^2=D^*D+\Ric+K^2\geq D^*D,
\end{align*}

we have that both these Riesz transforms are bounded in $L^2$ by 1.

Now we want to prove that, for a suitable range of values of $p$,
there exists a constant $C>0$ such that the inequality

\begin{align*}
\|\Hess \varphi\|_{L^p}\leq C\{\|\Delta \varphi\|_{L^p}+\|\varphi\|_{L^p}\}
\end{align*}

holds for every $\varphi\in C_c^{\infty}(M)$. Clearly, it is enough
to show that

\begin{align*}
\|Dd\varphi\|_{L^p}\leq C\|(-\Delta +K^2)\varphi\|_{L^p},\quad \text{on}\ \ C_c^{\infty}(M).
\end{align*}

which, in turn, is equivalent to

\begin{align*}
\|Dd(-\Delta +K^2)^{-1}\varphi\|_{L^p}\leq C\|\varphi\|_{L^p}.
\end{align*}

Using the spectral calculus we write

\begin{align*}
(-\Delta +K^2)^{-1}=(-\Delta+K^2)^{-\frac{1}{2}}(-\Delta+K^2)^{-\frac{1}{2}}
\end{align*}

On the other hand, it can be deduced from [?] that the following
commutation rule holds:

\begin{align*}
\end{align*}

Whence, we are reduced to prove that

\begin{align*}
\|D(\Delta_1+K^2)^{-\frac{1}{2}}d(-\Delta+K^2)^{\frac{1}{2}}\varphi\|_{L^p}
\leq C\|\varphi\|_{L^p}
\end{align*}

and this is implied by the existence of some constants $C_1, C_2>0$
such that

\begin{align*}
&\|d(-\Delta +K^2)^{-\frac{1}{2}}\varphi\|_{L^p}\leq C_1\|\varphi\|_{L^p};\\
&\|D(\Delta_1+K^2)^{-\frac{1}{2}}\omega\|_{L^p}\leq C_2\|\omega\|_{L^p}.
\end{align*}

We have the following estimates for balls $B$ of $M$ centered in
$4B_j$ and with radii less than 8, $x\in B\cap 4B_j$ and functions $f$
in $L^2(M)$ with support in $B_j$:

1) The map $f\rightarrow \SM^{\#}_{4B_j,\tilde{T},n}f$ is bounded from
$L^p(B_j,\mu)$ into $L^p(4B_j,\mu)$ for all $p>2$, where

\begin{align*}
\SM^{\#}_{4B_j,\tilde{T},n}f(x)=\dsup_{B \, ball\ in\  M,\, B\ni x}\lf(\frac{1}{\mu(B\cap 4B_j)}\dint_{B\cap 4B_j}
\lf|\wz{T}(I-P_{r^2})^nf(y)\r|^2\,\mu(dy)\r)^{1/2}.
\end{align*}

2) If $p>2$, then for $k\in \{1,2,\cdots,n\}$,

\begin{align}\label{eq2}
\left(\frac{1}{\mu(B\cap 4B_j)}\int_{B\cap 4B_j}|\wz{T}P_{kr^2}f|^p\,d\mu\right)^{1/p}\leq C(\SM_{4B_j}(|\wz{T}f|^2)+(S_jf)^2)^{1/2}(x)
\end{align}

\begin{align}\label{eqn-squref}
(S_jf)^2:=\SM_{\rm loc}(|\wz{T}f|^2\mathbbm{1}_{M\setminus 4B_j})+\SM_{\rm loc}(|\nabla h|^2)+\SM_{4B_j}(|f|^2).
\end{align}

§ APPENDIX

In this appendix, we include the proof
of (CZ)$(2)$ in the case that the underlying manifold $M$ has a lower bound Ricci curvature
for reader's convenience (see [26, 36]).
To be precise, let $M$ be a complete Riemannian manifold satisfying the curvature condition
(<ref>).  For any
$u\in C_0^{\infty}(M)$, Bochner's formula gives

\begin{align}\label{eqn-BI}
-\frac{1}{2}\Delta |\nabla u|^2=|\Hess u|_\HS^2-g(\nabla \Delta u, \nabla u)+\Ric(\nabla u,\nabla u).
\end{align}

Integration by parts, together with the lower Ricci bound, leads to

\begin{align*}
0&=-\frac{1}{2}\int_M \Delta |\nabla u|^2\\
&=\int_M |\Hess u|_\HS^2-\int_M g(\nabla \Delta u, \nabla u)+\int_M \Ric(\nabla u, \nabla u)\\
&\geq \int_M |\Hess u|_\HS^2-\int_M(\Delta u)^2-K\int_M |\nabla u|^2\\
&=\int_M |\Hess u|_\HS^2-\int_M(\Delta u)^2-K\int_M u\Delta u.
\end{align*}

Then, by estimating the last integral via Young's inequality, we
obtain for any $\epz>0$,

\begin{align}\label{T2-2}
\|\,|\Hess u|\,\|_{L^2}^2\leq\|\,|\Hess u|_\HS\,\|_{L^2}^2\leq \frac{K\epz^2}{2}\|u\|_{L^2}^2+\lf(1+\frac{K}{2\epz^2}\r)\|\Delta u\|_{L^2}^2,
\end{align}

which establishes CZ(2).

Inequality (<ref>) extends from $u\in C_0^{\infty}(M)$ to
$u\in H^{2,2}(M)$.  Thus, in particular, if $u\in L^2(M)$ is a
distributional solution to the Poisson equation

\begin{align}\label{eqn-Poisson}
\Delta u=f \quad \text{on }M
\end{align}

for some $f\in L^2(M)$, then (<ref>) provides an $L^2$-Hessian
estimate of the solution $u$.

Recall that for $L^2$-gradient estimates of $u$ in
(<ref>), the lower bound on the Ricci tensor is not
needed.  Indeed, integrating (<ref>) by parts and using
Young inequality, one obtains directly for every $\epz>0$,

\begin{align*}
\int_M|\nabla u|^2=\int_Mf u
\leq \frac{\epz^2}{2}\int_M u^2+\frac{1}{2\epz^2}\int_M f^2,
\end{align*}

which is the $L^2$-gradient estimate

\begin{align}\label{eqn-Lge}
\|\nabla u\|^2_{L^2}\leq \frac{\epz^2}{2}\|u\|_{L^2}^2+\frac{1}{2\epz^2}\|f\|_{L^2}^2.
\end{align}

[1]
Michael Aizenman and Barry Simon, Brownian motion and Harnack
inequality for Schrödinger operators, Comm. Pure Appl. Math.
35 (1982), no. 2, 209–273. 644024

[2]
Marc Arnaudon, Holger Plank, and Anton Thalmaier, A Bismut type formula
for the Hessian of heat semigroups, C. R. Math. Acad. Sci. Paris
336 (2003), no. 8, 661–666. 1988128

[3]
Pascal Auscher, Thierry Coulhon, Xuan Thinh Duong, and Steve Hofmann,
Riesz transform on manifolds and heat kernel regularity, Ann. Sci.
École Norm. Sup. (4) 37 (2004), no. 6, 911–957. 2119242

[4]
Robert Baumgarth, Baptiste Devyver, and Batu Güneysu, Estimates for
the covariant derivative of the heat semigroup on differential forms, and
covariant Riesz transforms, arXiv:2107.00311 (2021).

[5]
Jean-Michel Bismut, Large deviations and the Malliavin calculus,
Progress in Mathematics, vol. 45, Birkhäuser Boston, Inc., Boston, MA,
1984. 755001

[6]
Alberto P. Caldéron and Antoni S. Zygmund, On the existence of
certain singular integrals, Acta Math. 88 (1952), 85–139.

[7]
Gilles Carron, Riesz transform on manifolds with quadratic curvature
decay, Rev. Mat. Iberoam. 33 (2017), no. 3, 749–788. 3713030

[8]
Isaac Chavel, Eigenvalues in Riemannian geometry, Pure and Applied
Mathematics, vol. 115, Academic Press, Inc., Orlando, FL, 1984, Including a
chapter by Burton Randol, With an appendix by Jozef Dodziuk. 768584

[9]
Xin Chen, Xue-Mei Li, and Bo Wu, Logarithmic heat kernels: estimates
without curvature restrictions, arXiv:2106.02746 (2021).

[10]
Bennett Chow and Richard S. Hamilton, Constrained and linear Harnack
inequalities for parabolic equations, Invent. Math. 129 (1997),
no. 2, 213–238. 1465325

[11]
Ronald R. Coifman and Guido Weiss, Analyse harmonique non-commutative sur
certains espaces homogènes, Lecture Notes in Mathematics, Vol. 242,
Springer-Verlag, Berlin-New York, 1971, Étude de certaines intégrales
singulières. 0499948

[12]
Thierry Coulhon and Xuan Thinh Duong, Riesz transforms for $1\leq p\leq
2$, Trans. Amer. Math. Soc. 351 (1999), no. 3, 1151–1169.

[13]
, Riesz transforms for $p>2$, C. R. Acad. Sci. Paris Sér. I
Math. 332 (2001), no. 11, 975–980. 1838122

[14]
, Riesz transform and related inequalities on noncompact
Riemannian manifolds, Comm. Pure Appl. Math. 56 (2003), no. 12,
1728–1751. 2001444

[15]
E. Brian Davies, The state of the art for heat kernel bounds on
negatively curved manifolds, Bull. London Math. Soc. 25 (1993),
no. 3, 289–292. 1209255

[16]
, Non-Gaussian aspects of heat kernel behaviour, J. London
Math. Soc. (2) 55 (1997), no. 1, 105–125. 1423289

[17]
Xuan T. Duong and Derek W. Robinson, Semigroup kernels, Poisson bounds,
and holomorphic functional calculus, J. Funct. Anal. 142 (1996),
no. 1, 89–128. 1419418

[18]
K. David Elworthy and Xue-Mei Li, Formulae for the derivatives of heat
semigroups, J. Funct. Anal. 125 (1994), no. 1, 252–286.

[19]
, Bismut type formulae for differential forms, C. R. Acad. Sci.
Paris Sér. I Math. 327 (1998), no. 1, 87–92. 1650216

[20]
David Gilbarg and Neil S. Trudinger, Elliptic partial differential
equations of second order, Classics in Mathematics, Springer-Verlag, Berlin,
2001, Reprint of the 1998 edition.

[21]
Alexander Grigor’yan, Gaussian upper bounds for the heat kernel on
arbitrary manifolds, J. Differential Geom. 45 (1997), no. 1,
33–52. 1443330

[22]
Batu Güneysu, Nonrelativistic hydrogen type stability problems on
nonparabolic 3-manifolds, Ann. Henri Poincaré 13 (2012), no. 7,
1557–1573. 2982633

[23]
, Covariant Schrödinger semigroups on Riemannian
manifolds, Operator Theory: Advances and Applications, vol. 264,
Birkhäuser/Springer, Cham, 2017. 3751359

[24]
, Heat kernels in the context of Kato potentials on arbitrary
manifolds, Potential Anal. 46 (2017), no. 1, 119–134. 3595965

[25]
Batu Güneysu and Diego Pallara, Functions with bounded variation on a
class of Riemannian manifolds with Ricci curvature unbounded from below,
Math. Ann. 363 (2015), no. 3-4, 1307–1331. 3412360

[26]
Batu Güneysu and Stefano Pigola, The Calderón-Zygmund
inequality and Sobolev spaces on noncompact Riemannian manifolds, Adv.
Math. 281 (2015), 353–393. 3366843

[27]
, $L^p$-interpolation inequalities and global Sobolev
regularity results, Ann. Mat. Pura Appl. (4) 198 (2019), no. 1,
83–96, With an appendix by Ognjen Milatovic. 3918620

[28]
Batu Güneysu and Olaf Post, Path integrals and the essential
self-adjointness of differential operators on noncompact manifolds, Math. Z.
275 (2013), no. 1-2, 331–348. 3101810

[29]
Steve Hofmann and José María Martell, $L^p$ bounds for Riesz
transforms and square roots associated to second order elliptic operators,
Publ. Mat. 47 (2003), no. 2, 497–515. 2006497

[30]
Shouhei Honda, Luciano Mari, Michele Rimoldi, and Giona Veronelli,
Density and non-density of $C_c^\infty\hookrightarrow W^{k,p}$ on
complete manifolds with curvature bounds, Nonlinear Anal. 211
(2021), Paper No. 112429, 26. 4268755

[31]
Tosio Kato, Schrödinger operators with singular potentials, Israel
J. Math. 13 (1972), 135–148 (1973). 333833

[32]
Xue-Mei Li, Stochastic differential equations on noncompact manifolds,
University of Warwick, Thesis (1992).

[33]
Ludovico Marini and Giona Veronelli, The $L^p$-Calderón-Zygmund
inequality on non-compact manifolds of positive curvature, Ann. Global Anal.
Geom. 60 (2021), no. 2, 253–267. 4291611

[34]
Lei Ni, The Poisson equation and Hermitian-Einstein metrics on
holomorphic vector bundles over complete noncompact Kähler manifolds,
Indiana Univ. Math. J. 51 (2002), no. 3, 679–704. 1911050

[35]
Lei Ni, Yuguang Shi, and Luen-Fai Tam, Poisson equation,
Poincaré-Lelong equation and curvature decay on complete Kähler
manifolds, J. Differential Geom. 57 (2001), no. 2, 339–388.

[36]
Stefano Pigola, Global Calderón-Zygmund inequalities on complete
Riemannian manifolds, arXiv:2011.03220v1 (2020).

[37]
Christian Rose and Peter Stollmann, Manifolds with Ricci curvature in
the Kato class: Heat kernel bounds and applications, Analysis and Geometry
on Graphs and Manifolds, London Mathematical Society Lecture Note Series,
vol. 461, Cambridge University Press, 2020, pp. 76–94.

[38]
Barry Simon, Schrödinger semigroups, Bull. Amer. Math. Soc. (N.S.)
7 (1982), no. 3, 447–526. 670130

[39]
Daniel W. Stroock, An estimate on the Hessian of the heat kernel,
Itô's stochastic calculus and probability theory, Springer, Tokyo, 1996,
pp. 355–371. 1439536

[40]
Anton Thalmaier and Feng-Yu Wang, Derivative estimates of semigroups and
Riesz transforms on vector bundles, Potential Anal. 20 (2004),
no. 2, 105–123. 2032944

[41]
Caitlin Wang, The Calderón-Zygmund inequality on a compact
Riemannian manifold, Pacific J. Math. 217 (2004), no. 1,
181–200. 2105773

[42]
Feng-Yu Wang, Analysis for diffusion processes on Riemannian
manifolds, Advanced Series on Statistical Science & Applied Probability,
vol. 18, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2014.