Consider a collection of positive definite matrices $\{A_1,...,A_n\}\in\mathbb{S}_{++}^d$, the Karcher mean of these matrices is given by (see (5.4) in [1]):
$$
\min_{X\in\mathbb{S}_{++}^d} f(X):=\frac{1}{n}\sum_{i=1}^{n}\|\ln(X^{-1/2}A_i X^{-1/2})\|_F^2
$$
where $f$ is simply the sum of squared geodesic distances of $X$ to all the $A_i$s.

The Riemannian gradient of $f$ is given as (5.5) in [1]:
$$
\mathrm{grad} f(X) = \frac{1}{n}\sum_{i=1}^{n} X^{1/2}\ln(X^{1/2}A_i^{-1}X^{1/2})X^{1/2}
$$

My question is how to calculate the Riemannian Hessian of this function $f$, if it enjoys closed form? 

(A more generic question is that can the Riemannian Hessian of squared geodesic distance being expressed in closed form by operators on Riemannian manifolds, such as exponential and logarithm mappings?)

Any comment or reference are appreciated!

References:
[1] O. P. Ferreira, M. S. Louzeiro, and L.F. Prudente. "Gradient method for optimization on Riemannian manifolds with lower bounded curvature." SIAM Journal on Optimization 29.4 (2019): 2517-2541.