Let $M$ be a compact Riemann Surface. Let $P$ be a fixed point on $M$ and let $\delta_{P}$ be the Dirac point distribution at $M$. Consider the fundamental solution of the heat equation $$ (\partial_{t}-\Delta)K(x,y,t)=0, K(x,y,0)=\delta_{x}(y) $$ and its solution $$ f(x)=e^{\Delta}\delta_{P}(x)=\int K(x,y,1)\delta_{P}(y)dy=K(x,P,1) $$ I am wondering if there is any special meaning attached to the case when $t=1$ and $M$ is taken as a homogeneous space. The idea seems natural that one start Brownian motion with a pointed probability distribution and let time flows, then $f(x)$ can be used as a measure of the "distance" between $x$ and $P$. In the case $M=\mathbb{S}^{1}$ via its universal cover $\mathbb{R}^{1}$ we recovered a sum of Gaussians that resembles theta function.
However when I look up online, the explicit formulas we have are rather difficult to interpret. For $g\ge2$ they are quotients of $\mathbb{H}$, and thanks to Mckeans' work we know $$ K(x,y,t)=\frac{\sqrt{2}}{(4\pi t)^{3/2}}e^{-t/4 }\int^{\infty}_{p}\frac{se^{-s^2/4t}}{(\cosh(s)-\cosh(p))^{1/2}}ds,p=d_{\mathbb{H}}(x,y) $$ I do not really know any good interpretation of this when we fix $y=P,t=1$ and let $x$ vary then add up the copies corresponding to $\pi_{1}(M)$. In particular, I am curious what does this say related to the spectrum of the Laplacian in $M$. It should relate to analytic torsion of the Laplacian on the manifold, but so far I have not found any way to relate the two objects.
This coming from the following observation in Arakelov theory: In the arithemetic surface setting to obtain a Riemann-Roch formula for a line bundle $L$ is the same thing as giving a metric for the determinant line bundle $\lambda(L)=\lambda(H^{0}(L))\otimes \lambda(H^{1}(L)^{*})$. I am wondering if there is anyway to define $H^{1}$ directly which satisfies Serre duality. In one dimensional case, the regularized $h^{1}$ proposed by Van Der Geer and Rene Schoof is essentially a weighted sum of Gaussian functions, which may be interpreted as the probability density of standard Brownian motion on the real line covering the interval $[0,1]$. However in two dimensional situation it is not entirely self-evident that analgous methods give us a self-dual object. A standard choice in the community is to work with the Quillen metric, but it is not clear to me how this may relate to Quillen metric either. The question at here is really about how to understand this probability density better, as the original problem seems still totally out of reach.
There is nothing "special" about $t=1$. Indeed if we rescale $M$ or change $t=1$ with $t=s$, then the corresponding heat kernel would also change. For me the most direct reason may be in $\mathbb{R}^{1}$ case for $t=1$ the Fourier transform of the heat kernel is itself. I am not sure what the Fourier transform of the heat kernel on the hyperbolic space is for fixed $t$ with respect to $p$, for example. And in general this self-duality should stem from some form of Poincare duality, as Serre duality can be viewed as Poincare duality for Riemann surfaces.
The reason I viewed it as a sort of "distance" is again motivated by Arakelov theory. In the Vander Geer-Schoof paper they used an effectivity function of the form (see page 6) $$ e((f)+D)=e^{-\pi |f|^{2}_{D}} $$ So to mimic the behavior of effectivity function in $\mathbb{R}$ it seems natural to consider its counterpart on the Riemann surface, such that for points "closer" to $P$ the value $f(x)$ should be greater as it got hitted more often. Some numerical computation showed the function $f_{P}(x), P=0$ indeed decays very fast with respect to $p=d(x,0)$. However I am unable even to plot it with respect to $p$.