Sharp Gaussian upper bounds on Heat Kernel 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. 
 A: If $\mathbb{M}$ is a compact manifold, the bounds you are interested in are a consequence of the Minakshisundaram-Pleijel expansion
$p_t(x,y)\sim_{t \to 0} \frac{e^{-\frac{d^2(x,y)}{4t}}}{(4\pi t)^{n/2}}\sum_{k=0}^{+\infty} u_k(x,y)t^k $
which holds uniformly on compact sets of $\mathbb{M}\times \mathbb{M} -C_{\mathbb{M}}$ where 
$C_\mathbb{M}=\{ (x,y), x\text{ is in the cut-locus of }y\}$
A non-probabilistic proof of it relies on the parametrix method.
A classical reference is Chapter 6 in Chavel:Eigenvalues in Riemannian geometry. You will find the uniform asymptotic in small times page 154. Another reference is Chapter 1 in Berline-Getzler-Vergne: Heat kernels and Dirac operators. 
If the manifold is non-compact, an additional  cutoff argument is needed to obtain the uniform expansion. The quick one I know uses probability. It may be found in Chapter 5 of the book by Hsu: Stochastic analysis on manifolds.
Edit.
As you  noted, there is a mistake  in the proof by Chavel (interestingly this mistake seems to originate from Donnelly which is cited in Chavel)... As far as I can see now, the proof by Chavel can be fixed to only give the estimate
$\left| p_t(x,y)- \frac{e^{-\frac{d^2(x,y)}{4t}}}{(4\pi t)^{n/2}}\sum_{k=0}^{N} u_k(x,y)t^k \right| \le C_{K,N}t^{N+1-n/2}  e^{-\frac{d^2(x,y)}{(4+\varepsilon)t}}$
However we can prove by probability methods  that
$\left| p_t(x,y)- \frac{e^{-\frac{d^2(x,y)}{4t}}}{(4\pi t)^{n/2}}\sum_{k=0}^{N} u_k(x,y)t^k \right| \le C_{K,N}t^{N+1-n/2}  e^{-\frac{d^2(x,y)}{4t}}$
So you get the bounds you are interested in. Observe that no curvature conditions are needed.
If you want  references, the result by Molchanov was reproved and generalized by Robert Azencott in this paper, by Jean-Michel Bismut in this book , and in this paper by Gerard Ben Arous. A more recent reference is this paper. You can have a look to the references.
A last comment is that if you assume that Ricci is bounded from below by -K, then we have on any compact subset of $\mathbb{M}\times\mathbb{M}$ and $t \in [0,T]$,
$p_t(x,y)\ge \frac{C}{t^{n/2}} e^{-\frac{d^2(x,y)}{4t}}$
This is a consequence of the parabolic Harnack inequality (no probability here !). You can see the proof on my blog, see the last proposition.
