 Solution

Retain $\Delta =\mathrm{div}\nabla$ and define the heat operator by $L={\partial }_t-{\Delta }_q$. A Riemannian heat kernel $H(p,q,t)$ satisfies $L_{q}H=0$ for $t>0$ and has the delta initial limitOn standard Euclidean space the kernel isOn the Riemannian manifold, the Gaussian function factor in the local construction uses $r(p,q)=\mathrm{dist}(p,q)$ in place of $|p-q|$.

Choose an open neighbourhood $\mathcal{U}$ of the diagonal such that every pair $(p,q)\in \mathcal{U}$ is joined by a unique short minimizing geodesic and $q={\exp }_{p}x$ with the inverse exponential map smooth in both variables. It can be obtained from normal coordinates tubes of radius smaller than the local injectivity radius, or from normal neighbourhoods with geodesic convexity. Then $r^2$ is smooth on $\mathcal{U}$. Write the Riemannian volume form in normal coordinates as $d{V}_q=J(p,x)dx$, where $J>0$ is smooth and $J(p,0)=1$.

For fixed $p$, radial divergence in normal coordinates givesFor $g=(4\pi t{)}^{-d/2}{e}^{-r^2/(4t)}$, use $|\nabla r|=1$ away from the diagonal and the product rule to obtainThe identities extend smoothly in the expressions involving $r{\partial }_r$, even though $r$ itself is not smooth on the diagonal. The singular powers cancel ifThe correct initial normalization is $w_0(p,p)=1$. These heat-kernel transport equations have the explicit solutionsand, recursively,To derive the integral, multiply the radial equation by $r^{i-1}{J}^{1/2}$ and integrate the derivative of $r^i{J}^{1/2}{w}_i$ from zero to $r$. The zero integration constant is forced by smoothness when $i\ge 1$. Changing variables from radial distance to $s\in [0,1]$ gives the displayed nonsingular expression. Its integrand is smooth in $(p,x,s)$, so differentiation under the integral proves inductively that every $w_i$ is smooth across the diagonal. In particular $w_i(p,p)=({\Delta }_q{w}_{i-1})(p,p)/i$.

Summing the product-rule identity makes adjacent terms telescope. For $S_k=g{\sum }_{i=0}^k{t}^i{w}_i$, the only term remaining isThis also fixes the sign convention: with the nonnegative Laplace-Beltrami operator $P=-\Delta$, the right-hand side would instead be $+g{t}^k{P}_q{w}_k$.

A heat parametrix is an approximate kernel with the same delta initial limit and an error under $L$ that is sufficiently regular, or flat, at $t=0$ to be corrected by time convolution. On a closed compact manifold choose a smooth cutoff $\chi (p,q)$ equal to one near the diagonal and supported in $\mathcal{U}$. A finite-order explicit parametrix isNear the diagonal its error is $-\chi g{t}^N{\Delta }_q{w}_N$. Derivatives of $\chi$ produce errors supported a positive distance from the diagonal, where $e^{-r^2/(4t)}$ decays faster than every power of $t$. Taking $N$ large makes the whole error as regular at time zero as any specified finite number of derivatives requires. Its initial delta limit follows from $J(p,0)=w_0(p,p)=1$ and the Euclidean Gaussian function change of variables $x=\sqrt{t}y$.

One may obtain a smooth error flat at zero by summing all coefficients with time cutoffs. Let $\rho \in C^{\infty }([0,\infty ))$ equal one near zero and zero for arguments at least one, and choose ${ϵ}_i\downarrow 0$ sufficiently fast. SetFor each $t>0$ the sum is locally finite. Choosing ${ϵ}_i$ successively so that the $i$th term and its derivatives through order $\lfloor i/2\rfloor$ have bounds $2^{-i}$ on their cutoff transition regions gives the usual smooth asymptotic sum. For any fixed derivative order, its tail is then smaller than an arbitrarily high power of $t$; the transport identities cancel all earlier orders. Consequently $R=L_{q}P$ extends smoothly and flatly to $t=0$.

Here is the Volterra parametrix correction, with the order of the Volterra convolution of kernels factors specified. DefineSince $L$ acts in the second spatial variable and $P$ has the delta initial limit, $L_q(Q\ast P)=Q+Q\ast R$. ThusIndeed, $Q+R+Q\ast R=0$ by the geometric-series cancellation, whence $L_{q}H=0$. If $|R|\le C$ on a bounded time interval and $v=\mathrm{vol}(M)$, thenThis proves convergence; the smooth flat error permits the corresponding derivative estimates. The correction has zero initial limit, so $H$ retains the delta initial condition. Uniqueness of the heat equation on a closed compact manifold identifies this kernel with the global Riemannian heat kernel.

The local coefficient construction requires no compactness. The global uniformly bounded convolution argument just given uses a closed compact manifold, the usual spectral-geometry setting. For noncompact or incomplete Riemannian manifolds one must specify the heat realization and justify global convergence separately; a canonical choice is the minimal kernel obtained as the increasing limit of Dirichlet kernels on a smooth relatively compact exhaustion. This gives the minimal heat semigroup rather than asserting an unqualified uniqueness statement at infinity or at a missing boundary.