Hodge laplacian of distance function Let $p$ be a given point on a Riemaniann manifold $\mathcal{M}$. The distance function to point $p$ is denoted $f_p$ :
$$ f_p(q) = \operatorname{dist}(p,q)$$
The exterior derivative is denoted $\mathrm{d}$ and the codifferential is denoted $\delta$. Then the Hodge laplacian is $\Delta = \mathrm{d}\delta + \delta \mathrm{d} = (\mathrm{d} + \delta)^2$ which reduces to $\delta \mathrm{d}$ in the case of functions.
My questions are the following :


*

*I think it can be shown that $f_p^2$ (the squared distance function to point $p$) is $C^\infty$ in a neighborhood of $p$. Can anyone confirm ?

*What is the value of the Hodge laplacian of $f_p^2$ ? By analogy with the case $\mathcal{M} = \mathbb{R}^n$, is my guess that it is a constant function right ? Do we have $\Delta f_p^2 = 2n$  where $n$ is the dimension of $\mathcal{M}$ ?

 A: 1) By first removing the cut locus of $p$, the remaining subset of $M$ is an open subset diffeomorphic to some open subset of $T_p M$. Furthermore, if you place yourself in spherical normal coordinates $(r, \sigma)$ (where $\sigma$ encodes all the angular variables), the function $d_p ^2$ becomes $r^2$, which is obviously smoooth.
For a rigorous statement (but following the same idea), see theorem 3.6 at page 166, in chapter IV of "Foundations of Differential Geometry", volume 1, by Kobayashi and Nomizu.
2) In the same spherical normal coordinates, the Laplacian looks like
$$\Delta f \ (q) = \frac {\partial ^2 f} {\partial r^2} (q) + H(p,q) \frac {\partial f} {\partial r} (q) + \frac 1 {r^2} \Delta_S f \ (q) ,$$
where $r = d(p,q)$, $S$ is the geodesic sphere of radius $r$ centered at $p$, $\Delta_S$ is the Laplacian restricted to $S$ and $H(p,q)$ is the mean curvature at $q$ of $S$. It follows that
$$\Delta d_p ^2 \ (q) = \Delta r^2 = 2 + H(p,q) 2r = 2 + 2 H(p,q) d_p (q) .$$
In $\Bbb R^n$, the mean curvature of the sphere of radius $r$ is $H = \frac {n-1} r$, so the above becomes $2n$. Otherwise, as your intuition hinted at, the curvature does play a role and makes the result non-constant, in general.
For an alternative formulation of the Laplacian see also page 48 of "Spectral Theory and Geometry" edited by E.B. Davies and Y. Safarov.
A: The answer to your first question is yes.  You can view the distance function you defined 
$$f_p(q) = d(p, q)$$
as the radial component of the inverse of the exponential map at $p$.  The exponential map is known to be a diffeomorphism in a neighborhood $p$.  You can define the full neighborhood in which it is smooth in terms of its conjugate points according to this Mathoverflow Thread:

The distance function is differentiable at $(p,q) \in M \times M$ if and only if there is a unique length-minimizing geodesic from p to q. Furthermore, the distance function is $C^\infty$ in a neighborhood of $(p,q)$ if and only if $p$ and $q$ are not conjugate points along this minimizing geodesic.

Regarding the answer to your second question, I can't find a specific counter-example but I see no reason why it would be true.  I did do various searches of literature and found references to bounding the Laplacian of distance functions.  e.g. In On the differential structure of metric measure spaces and applications it is stated on page $2$ that 
$$\Delta f_p(q) \leq \frac{n - 1}{f_p(q)}$$ 
when $M$ is a manifold with non-negative Ricci curvature.
