I am looking for references (with proof) for the following statement:

Let $(M, g)$ be a Riemannian manifold with bounded curvature and let $p_t(x , y)$ be the heat kernel of $M$. Let $K$ be compact in $M \times M$ such that $K$ is disjoint from the cut locus and let $T>0$. Then there exist constants $c, C > 0 $ (depending on $(M, g)$, $K$ and $T$) such that $$\frac{c}{t^{n/2}} e^{-\frac{1}{4t}d(x, y)^2} \leq p_t(x, y) \leq \frac{C}{t^{n/2}} e^{-\frac{1}{4t}d(x, y)^2}$$ uniformly for all $t \in (0, T]$ and $(x,y) \in K$.

There are many variants of this theorem in the literature: Either with a different power of $t$ in the denominator on the right-hand-side, or with a denominator greater than $4t$ in the exponent, but valid for all points in $M\times M$ (also for cutpoints) or for all times (papers and books by Grigor'yan, Davies and many others).

The only source for this exact statement (sort of) I could find was the paper "Diffusion processes and Riemannian Geometry" by Molchanov, where the proof is based on stochastic analysis (which is hard for me to understand), and also he just writes $$p_t(x, y) \sim \frac{J(x, y)^{-1/2}}{4\pi t^{n/2}} e^{-\frac{1}{4t}d(x, y)^2}$$ ($J(x, y)$ is the Jacobian determinant of the Riemannian exponential map), but never explains what the symbol $\sim$ actually means. And since I do not understand the proof properly, I am not even sure that it is supposed to mean what I wrote above.

Also it would be nice (for me) to have a proof that does not rely on stochastic analysis.

/Edit: This is an answer to Fabrice Boudoin's Post which is too long for a comment. The problem is the following: Berline-Getzler-Vergne and Berger-Gauduchon-Mazét proof that the asymptotic heat kernel expansion satisfies
$$ \Bigl| p_t(x, y) - (4\pi t)^{-n/2}e^{-\frac{1}{4t}d(x, y)^2}\sum_{j=0}^N t^j \Phi_j(x, y)\Bigr| \leq C_{K, N}t^{-n/2 +N+1}$$
for all $N$, uniformly on compact sets $K$ disjoint from the cut locus. However, this does *not* imply the statement that I am looking for.

The statement that I mean *would* be implied by the statement that
$$ \Bigl| \frac{p_t(x, y)}{(4\pi t)^{-n/2}e^{-\frac{1}{4t}d(x, y)^2}} -\sum_{j=0}^N t^j \Phi_j(x, y)\Bigr| \leq C_{K, N}t^{N+1}.$$
Chavel has this statement, but his proof relies on Lemma 1 on page 152, which is wrong (and he does not include a proof anyways). For the proof, it is refered to Berger-Gauduchon-Mazet, who do not prove Lemma 1 (which is not surprising as it is wrong...). The Book of Hsu you mentioned does contain the statement, but he just refers to Chavel for a proof.

compactneighbourhood of diagonal exist for non-compact $M$? – R W May 29 '14 at 2:04