Normal-coordinate divergence-form heat parametrix ($B_k(x)=-{\int }_0^1{s}^{k-1}(L{B}_{k-1})(sx)ds$)

For $L=-{\partial }_j(g^{ij}{\partial }_i)$ in geodesic normal coordinates, the Gauss lemma gives $g^{ij}{x}_i=x_j$. Substituting $F=(4\pi t{)}^{-n/2}{e}^{-|x{|}^2/(4t)}{\sum }_{k\ge 0}{t}^k{B}_k$ into $({\partial }_t+L)F=0$ therefore yields $(x\cdot \nabla +k)B_k=-L{B}_{k-1}$, with $B_0=1$. Integration along radial rays gives the displayed formula. Every coefficient is smooth at the centre. Any smooth solution of $(x\cdot \nabla +k)h=0$ is zero for $k>0$, since $s^{k}h(sx)$ is constant and tends to zero at $s=0$. This proves uniqueness. In this special divergence-form case $L1=0$, so every $B_k$ for $k>0$ is zero: the Gaussian itself solves the equation exactly on the coordinate ball. The positive Laplace-Beltrami operator instead has a metric volume-density term and a nonconstant leading amplitude.