For the integral kernel of the Laplacian $\Delta$ on $\mathbb{R}^n$, consider the resolvent $R(\lambda) := (\lambda - \Delta)^{-1}$ and let $R(\lambda; x, y)$ be its kernel, which is a smooth function away from the diagonal $\{x=y\}$. If I calculated correctly, we have estimates of the form $$|R(\lambda;x, y)| \leq C e^{-\delta(\lambda)|x-y|},$$ whenever $|x-y| \geq 1$, where $\delta(\lambda)$ is a certain positive constant, depending (quite explicitly) on the distance of $\lambda$ to the spectrum of $\Delta$.
Q: On general Riemannian manifolds, do we still have estimates like this on the resolvent kernel of the Laplacian (or more general operators)?
More specifically, if $M$ is a complete Riemannian manifold, say with bounded geometry, and $\Delta$ is the Laplace-Beltrami operator, is it true that for any $\varepsilon>0$, there exist constants $C, \delta>0$ such that we have $$|R(\lambda;x, y)| \leq C e^{-\delta\,d(x, y)}$$ whenever $d(x, y) \geq 1$ and $\lambda$ has distance at least $\varepsilon$ from the spectrum of $\Delta$? Here $d(x, y)$ is the Riemannian distance function.