I am trying to gather some intuition on the connection on diffusion processes and Riemannian geometry, with only very limited knowledge of the latter.
First, let us consider a Brownian Motion in euclidean space. For $\Sigma$ a positive definite matrix, $W_t$ standard $n-$dimensional brownian motion, $X_t$ solving the stochastic differential equation $$ d X_t = \Sigma^{\frac{1}{2}} d W_t$$
we have that the density $p(x,y,t)$ for $X_t$ starting at $y$ solves the Fokker-Planck equation $$ \frac{d}{dt} p = \frac{1}{2} \Sigma_{ij} \frac{\partial}{\partial x_i} \frac{\partial}{\partial x_j} p$$ where we use Einstein summing convention.
The density is given by $p(x,y,t) \propto \exp \left(\frac{-(x-y)^T \Sigma^{-1}(x-y)^T} {2t}\right)$ (where the proportionality is since we still require normalization).
Now some geometry: Let $(M,g)$ be a Riemannian manifold (we assume compactness and everything else we might need for the definitions to work).
The gradient $\operatorname{grad}$ is defined implicitly by obeying the relation $\langle \operatorname{grad} f, X \rangle_g = d f(X)$ for any vector vield $X$. The adjoint of the gradient is the divergence $\operatorname{div}$. One then has the Laplace-Beltrami-operator on the manifold given by $$\Delta_M = \operatorname{div}\operatorname{grad}.$$
One can explicitly write the Laplace-Beltrami operator as $$\Delta_M f = \frac{1}{\sqrt{\det g}} \frac{\partial}{\partial x_j} (\sqrt{\det g} g^{ij} \frac{\partial}{\partial x_i} f)$$ where $g^{ij}$ refers to the inverse of the metric tensor.
We now see that our Fokker-Planck equation includes a special case of the Laplace-Beltrami-Operator, where the manifold is just euclidean space and the metric is the constant $\Sigma^{-1}$.
Indeed, in these notes the definition of a Brownian motion on a manifold is taken to be the process which has density $p(x,y,t)$, where $p$ is the minimal, positive solution to $$\frac{d}{dt}p = \frac{1}{2} \Delta_M p, \ p(x,y,0) = \delta_y(x).$$
Now on our manifold, we have the notion of a geodesic distance, which is given by
$$ d_M(x,y) = \inf \{ L(\gamma) | \gamma \text{ a piecewise continuos curve from } x \text{ to } y\}.$$
In our standard euclidean case this is just a straight line. In particular, the geodesic distance on $E = (\mathbb{R}^n,\Sigma^{-1})$ should be given by $$d_E(x,y) = \sqrt{(x-y)^T \Sigma^{-1} (x-y)}.$$
Now by comparison with our density for the Brownian motion, we see that $$\ln p(x,y,t) \propto - d_E(x,y)^2.$$
This makes of course intuitively sense to be the generalization of the isotropic case: Diffusion is characterized by equiprobable motion in all directions, so by introducing a different covariane matrix of diffusion tensor we assume that we are in a space with some different metric, but still the transition probability only depends on this metric.
My question is now, whether this property generalizes to the definition of Brownian motion on arbitrary Riemannian manifolds above. That is, do we have that $$\ln p(x,y,t) \propto - d_M(x,y)^2$$ on a general Riemannian manifold $M$?
Intuitively, this should hold. Unfortunately, Im not familiar enough with differential geometry to prove this.
What might be useful is that geodesics should fulfill the 'geodesic equation' in local coordinates
$$\frac{\partial x_k} {\partial t^2} = -\Gamma^{k}_{ij} \frac{\partial x_i}{\partial t}\frac{\partial x_j}{\partial t} $$ with $\Gamma^{k}_{ij}$ the Christoffel symbol of the second kind. This also somewhat appears in the equation for the Laplace-Beltrami operator since
$$ \Gamma^i_{ij} = \frac{1}{\sqrt{\det g}} \frac{\partial}{\partial x_j} (\sqrt{\det g})$$
but I don't know if this helps.
I would appreciate any comments on the question, also if you know some good references to further study these connections.