Assume $E$ is a connected riemannian manifold with geodesic metric space structure given by $d$ and $P$ is a probability measure over $E$ with Borel sigma-algebra given by this metric structure. Also assume that $P$ is associated with some random variable $\xi$.

Traditional concept of expectation is not valid anymore as expression $$ \mathbb{E}\xi = \int_E x \, \mathrm{d}P(x) $$ doesn't make sense. However, we can define $f$-barycenter of measure $P$ to be any minimizer of $$ \sigma(P,f) = \min_{p \in E} \int_E f(d(p,x)) \, \mathrm{d}P(x). $$ where $f$ is a non-negative increasing Borel-measurable function (which is totally not important for the question so just assume $f = \mathrm{id}$ or $f = \mathrm{id}^2$). Of course, if $E$ is not Hadamard the barycenter is not unique, so denote by $\mathrm{Bar}(P,f)$ set of all barycentres. If $\sigma(P,f) = \infty$ then $\mathrm{Bar}(P,f) = E$, so assume this is not the case.

Number $\sigma(P,f)$ bears some resemblance to conventional variance. However, I'am not satisfied with this as in case of finite dimensional vector space variance turns into covarians matrix and in more general cases into covariance operator. This means that this $\sigma(P,f)$ is just a 'statistical mechanics style' variance. So my question is:

What is appropriate analogue of covariance operator for random variables taking values in riemannian manifold?

My thoughts so far:

As classical covariance matrix measure expected deviation of random variable into given 'direction', so appropriate covariance $\Sigma^\xi$ for any $\mu \in \mathrm {Bar}(P,f)$ must act as 2-tensor on vectors in tangent space $T_\mu E$ (maybe I actualy should use cotangent spaces or mix). Using physical notation for tangent vector let $v_\mu = (\mu,v,(U,\mathtt{x})),w_\mu=(\mu,w,(U,\mathtt{x}))$ where $(U,\mathtt{x})$ is a chart such that $\mathtt{x}(\mu) = 0$. It must be possible to define $$ \Sigma^\xi_\mu(v_\mu,w_\mu) = K_{\mathtt{x}(\xi_{|U})}(v,w) $$ where $K_{x(\xi_{|U})}$ denotes covariance operator for random vector $ \mathtt{x}(\xi_{|U})$ and $\xi_{|U}$ is constraint of $\xi$ to $U$ with probability $\mathbb{P}(\xi_{|U} \in A ) = \frac{P(A \cap U)}{P(U)}$ for any Borel set $A$. I don't like this construction becouse it is local and does not convey all possible information.

Much better definition must be $$ \Sigma^\xi_\mu(v_\mu,w_\mu) = \| v\|\|w\|\int_E \arg \min_{t \in I_{\mu,v}} d(x,\gamma_{\mu,v}(t)) \arg \min_{t \in I_{\mu,w}} d(x,\gamma_{ \mu, w}(t)) \, \mathrm{d}P(x) $$

where $\gamma_{\mu,v} : I_{\mu,v} \to E$ is a meridian-like constant speed curve such that $\gamma_{\mu,v}(0) = \mu$ and $$ \gamma'_{\mu,v}(t) = \frac{v_{\gamma_{\mu,v}(t)}}{\| v_{\gamma_{\mu,v}(t)} \|} $$ with $I_{\mu, v} = (-a,b)$ having $a,b \in \mathbb{R}_+ \cup \{\infty\}$ to be maximal numbers such that definition make sense. It also needs to be specified that in case of looping curve or multiple geodesics we must take $t$ with minimal absolute value and I do not know how to choose sign if $d(x,\gamma_{\mu,v}(t)) = d(x,\gamma_{\mu,v}(-t))$. I think such curves must be called maximal geodesics.

If there is a way to fix and improve this definition? How would you define curves $\gamma_{\mu, v} ?$

I like the second definition much more as it seems to be global and retains spirit of barycentric theory. It is also unclear if there is any sense in extending $\Sigma^\xi$ to arbitrary points $p \not \in \mathrm{Bar}(P,f)$ producing covariance tensor field. Can it be the case that $\mathrm{Bar}(P,f)$ is always a subriamanian manifold of $E$ and $\Sigma^\xi$ has a special meaning for it as a tensor field.

Is there any sense in covariance tensor field?

What realy bags me is that if we make set of all well-behaved probability measures $\mathcal{P}(E)$ over $E$ into riemannian space $(\mathcal{P}(E),w)$ (with Wasserstein geometry for example) then according to Pistones's nonparametric information geometry tangent spaces $T_{P}(\mathcal{P}(E),w)$ must be Orlich spaces of random vectors satisfying certain conditions in case $E = \mathbb{R}^n$. On the other hand, I am struggling to grasp what will be tangent vectors $X \in T_{P}(\mathcal{P}(E),w) $, maximal geodesics $\gamma_{P, X}$ and covariance tensor $\Sigma^\Xi_\mu$ for a random manifold measure $\Xi$ admitting a barycenter $\mu$.

Can you give me any clues ?

I understand that this question grew to a very big size, so I will be glad to see any answers.