Given a domain $\Omega$ of a complete Riemannian manifold $\mathcal{M}$ and define $\mathcal{A}$ to be the Laplacian with Neumann boundary condition on $\Omega$. We prove that, under appropriate conditions, the corresponding heat kernel satisfies the Gaussian upper bound$$h(t,x,y)\leq \frac{C}{\left[V_\Omega(x,\sqrt{t})V_\Omega (y,\sqrt{t})\right]^{1/2}}\left( 1+\frac{d^2(x,y)}{4t}\right)^{\delta}e^{-\frac{d^2(x,y)}{4t}},\;\; t>0,\; x,y\in \Omega .$$Here $d$ is the geodesic distance on $\mathcal{M}$, $V_\Omega (x,r)$ is the Riemannian volume of $B(x,r)\cap \Omega$, where $B(x,r)$ is the geodesic ball of center $x$ and radius $r$, and $\delta$ is a constant related to the doubling property of $\Omega$.As a consequence we obtain analyticity...