You just need to check that $$ \int_M H_k(x, y) \mathrm{d}y = 1 + O(t^{k+1})$$ for all $x \in M$ and for all $k$. For example, this follows directly from the method of stationary phase (just take a geodesic chart around $x$ that is so large that the support fo $\eta$ is in its domain. Translate to $\mathbb{R}^n$ and use the method of stationary phase there.

Therefore, your first equation in fact reduces to $$(4 \pi t)^{-n/2} \int_0^tC \tau^{k-n/2}(1 + O(\tau^{k+1}))\mathrm{d} \tau \leq C_1 t^{k+1}$$ which is obviously true.

So you just estimated to roughly when you took the sup-norm. Estimating the $L^1$-norm does the job.

You just need to check that $$ \int_M H_k(x, y) \mathrm{d}y = 1 + O(t^{k+1})$$ for all $x \in M$ and for all $k$. For example, this follows directly from the method of stationary phase (just take a geodesic chart around $x$ that is so large that the support fo $\eta$ is in its domain. Translate to $\mathbb{R}^n$ and use the method of stationary phase there.

Therefore, your first equation in fact reduces to $$(4 \pi t)^{n/2} \int_0^tC \tau^{k-n/2}(1 + O(\tau^{k+1}))\mathrm{d} \tau \leq C_1 t^{k+1}$$ which is obviously true.

So you just estimated to roughly when you took the sup-norm. Estimating the $L^1$-norm does the job.