The idea to construct a heat kernel is first construct a parametric in a small neighbourhood. Then use a bump function to extend it. And do convolution iteratively. (Reference: Laplacian on a Riemannian manifold [ROSENBERG].)
My question is will it spread all over the whole manifold? It seems that the bump function cut it off. But this should not happen for a heat kernel that it just vanish outside a neighbourhood. Because physically, the heat should be able to conduct to everywhere.
This bump function is not supported on an arbitrary neighborhood, but one very specific to the construction of $H(t,x,y)$. Rosenberg's argument proceeds by constructing the parametrix $H_k(t,x,y)$ in a neighborhood of the diagonal $M_{\text{diag}}\subset M\times M$. Namely the neighborhood $U_\epsilon=\{(x,y)\in M\times M:d(x,y)<\epsilon\}$.
The idea is that the heat kernel $H(t,x,y)$ vanishes for small $t$, provided that $x\neq y$. In fact, you should expect it vanish like $t^{-m/2}\text{exp}\left(-\frac{d(x,y)^2}{4t}\right)$, when comparing it with the Euclidean heat kernel. For this reason, you don't expect points in $(0,\infty)_t\times (M\times M\setminus U_\epsilon)$ to contribute much to the value of the Heat kernel, and the parametrix construction relies on this to some degree. The parametrix $H_k$ does not agree pointwise with the true heat kernel $H$, but instead differs by some small error (see theorem 3.22)
What he proves is that $H(t,x,y)=H_k(t,x,y)+Q_k\ast H_k(t,x,y)$, where the error function $Q_k$ satisfies $|Q_k|<Ct^{k-m/2}$, and that $Q_k\in C^l$ as long as $k>l+\frac{n}{2}$. In particular, on a closed manifold, convolution with $Q_k$ defines a compact operator (it is in this sense that the error is "small").
