Given a domain
${\rm\Omega}$ of a complete Riemannian manifold
${\mathcal{M}}$, define
${\mathcal{A}}$ to be the Laplacian with Neumann boundary condition on
${\rm\Omega}$. We prove that, under appropriate conditions, the corresponding heat kernel satisfies the Gaussian upper bound
$$\begin{eqnarray}h(t,x,y)\leq \frac{C}{[V_{{\rm\Omega}}(x,\sqrt{t})V_{{\rm\Omega}}(y,\sqrt{t})]^{1/2}}\biggl(1+\frac{d^{2}(x,y)}{4t}\biggr)^{{\it\delta}}e^{-d^{2}(x,y)/4t}\quad \text{for}~t>0,~x,y\in {\rm\Omega}.\end{eqnarray}$$ Here
$d$ is the geodesic distance on
${\mathcal{M}}$,
$V_{{\rm\Omega}}(x,r)$ is the Riemannian volume of
$B(x,r)\cap {\rm\Omega}$, where
$B(x,r)$ is the geodesic ball of centre
$x$ and radius
$r$, and
${\it\delta}$ is a constant related to the doubling property of
${\rm\Omega}$. As a consequence we obtain analyticity of the semigroup
$e^{-t{\mathcal{A}}}$ on
$L^{p}({\rm\Omega})$ for all
$p\in [1,\infty )$ as well as a spectral multiplier result.