Covariant derivative of inverse exponential vector field Let $M$ be a Riemannian manifold and denote by $\exp_p(v)$ the exponential map at $p \in M$ applied to $v \in T_p M$. Let $q \in M$ be fixed and let $U \subset M$ be a neighborhood of $q$ such that, for each $p \in U$, the map $\exp_p$ is a diffeomorphism in a neighborhood of $exp_p^{-1}(q)$.
Now, the map $V: p \mapsto \exp_p^{-1}(q)$ is a well-defined vector field on $U$ ($q$ is fixed!). Is there a "nice" expression for the covariant derivative $\nabla_X V$ along a given vector field $X$?
What I have tried so far
I understand that $V$ is a vector field such that $V(p)$ points in the (geodesic) direction of $q$ and such that $|V(p)| = d(p,q)$. On flat $\mathbb{R}^n$, we have $V(p) = q-p$ and the covariant derivative is simply $\nabla_X V = -X$.
In the Riemannian case, we have at least $\langle \nabla_X V, V \rangle = \langle -X, V \rangle$ for any vector field $X$. Proof: Consider $c(s,t) = \exp_q(s \exp_q^{-1}(\gamma(t))) = \exp_{\gamma(t)}((1-s) \exp_{\gamma(t)}^{-1}(q))$. Then $\frac{\partial c}{\partial s}(1,t) = -V(\gamma(t))$ and $|\frac{\partial c}{\partial s}(s,t)|$ does not depend on $s$. Therefore,
$$
  |V(\gamma(t))|^2 = \int_0^1 |\frac{\partial c}{\partial s}(s,t)|^2\,ds
$$
and
$$
\begin{aligned}
\langle \nabla_{\gamma(t)} V, V \rangle
&= \frac{1}{2}\frac{d}{dt}|V(\gamma(t))|^2
= \frac{1}{2}\frac{d}{dt}\int_0^1 |\frac{\partial c}{\partial s}(s,t)|^2\,ds \\
&= \langle \frac{\partial c}{\partial t}(1,t), \frac{\partial c}{\partial s}(1,t) \rangle
= -\langle \dot\gamma(t), V \rangle,
\end{aligned}
$$
where we used that $\frac{\partial c}{\partial t}(0,t) = 0$, $\nabla_{\partial_s}\frac{\partial c}{\partial s}(s,t) = 0$ and $\nabla_{\partial_s}\frac{\partial c}{\partial t}(s,t) = \nabla_{\partial_t}\frac{\partial s}{\partial t}(s,t)$.
My hypothesis
From explicit calculations on the sphere, I found that $V(p) = -F(p) \nabla F(p)$ where $F(p) = d(p,q)$. Then
$$
\nabla_X V = -\langle \nabla F, X\rangle \nabla F - F \nabla_X(\nabla F).
$$
From $|\nabla F(p)|^2 = 1$ ($F$ is 1-Lipschitz by triangle inequality), it's clear that $\nabla_X(\nabla F)$ is orthogonal to $\nabla F$ - just as expected. Actually, this formula could probably be more explicit and I don't have a proof for the relation between $F$ and $V$ in the general case.
 A: The vector field $V$ is the gradient of the function $f = -\frac12 r^2$, where $r = d(q,\cdot)$, the distance from $q$. To see this, note that the definition $V_p = \exp_p^{-1}(q)$ implies $\exp_p(V_p) = q$ for each $p\in U$, which is to say that the curve $\gamma(t) = \exp_p(tV_p)$ is a geodesic that goes from $q$ to $p$ in time $1$, and its velocity at $p$ is $V_p$. Thus the reverse geodesic $\sigma(t) = \gamma(1-t)$ is a radial geodesic that goes from $q$ to $p$ in time $1$. Its velocity at time $1$ is $\sigma'(1) = -\gamma'(0) = -V_p$. Every unit-speed radial geodesic starting at $q$ has velocity equal to $\operatorname{grad} r$ by the Gauss lemma. But the geodesic $\sigma$ is not unit-speed; instead, it traverses a distance $r(p)$ in time $1$, so its speed is $r(p)$. The upshot is $$
\newcommand{\grad}{\operatorname{grad}}
V_p = - r(p) \grad r|_p = \operatorname{grad} (-\tfrac12 r^2)|_p.
$$
The $(1,1)$-tensor field $\nabla V$ is thus the $(1,1)$-Hessian of $-\tfrac12 r^2$; let's denote it by $\mathscr H_f = \nabla (\grad f)$. We can also view it as an endomorphism field $\mathscr H_f\colon \Gamma(TM) \to \Gamma(TM)$, which acts on a vector field $X$ by
$$
\mathscr H_f(X) = \nabla _X V = \nabla_X(\grad f).
$$
Another way to look at it is that $\mathscr H_f = \nabla(\operatorname{grad} f)$ is the ordinary Hessian $\nabla^2 f$ with one index raised. We can write it as
\begin{align*}
\mathscr H_f  &= \nabla\operatorname{grad}(-\tfrac12 r^2)\\
& = 
\nabla(-r \operatorname{grad}r) =
-\operatorname{grad} r \otimes dr - r \nabla \operatorname{grad} r\\
&= -\operatorname{grad} r \otimes dr - r\mathscr H_r.\tag{$*$}
\end{align*}
To understand its action on an arbitrary vector field $X$, it is convenient to decompose $X$ as $X = a \operatorname{grad} r + Y$, where $Y$ is orthogonal to $\grad r$ (or equivalently $dr(Y)=0$). Because the integral curves of $\grad r$ are geodesics, $\nabla_{\grad r}(\grad r)\equiv 0$, which implies that $\mathscr H_r$ annihilates $\grad r$. Thus the action of $\mathscr H_f$ on $\grad r$ is just the action of the $-\grad r\otimes dr$ term, which gives $\mathscr H_f(\grad r)=- \grad r$. So to understand the behavior of $\mathscr H_f$, it suffices to consider its action on vectors orthogonal to $\grad r$, and for this we just need to understand the action of $\mathscr H_r$.
There are explicit formulas for $\mathscr H_r$ in constant-curvature manifolds. For a metric with constant sectional curvature $c$ and a vector $Y\perp \grad r$,
$$
\mathscr H_r(Y) = t_c(r) Y,
$$
where 
$$
t_c(r) = 
\begin{cases}
\dfrac 1 r, & c = 0,\\
\dfrac 1 R \cot \dfrac r R, & c= \dfrac 1{R^2}>0,\\
\dfrac 1 R \coth \dfrac r R, & c= -\dfrac 1{R^2}<0.
\end{cases}
$$
Finally, if the sectional curvatures of $M$ are bounded above or below by a constant, then there are comparison theorems that relate $\mathscr H_r$ to its constant-curvature counterpart. If $M$ has sectional curvatures bounded above by a constant $c$, then $\langle\mathscr H_r(Y),Y\rangle \ge t_c(r)|Y|^2$ for $Y\perp\grad r$ in any normal neighborhood of $q$, and if the sectional curvatures are bounded below by $c$, then $\langle\mathscr H_r(Y),Y\rangle\le t_c(r)|Y|^2$ for $Y\perp\grad r$ on the set where $r<\pi R$. This is proved by observing that $\mathscr H_r$ satisfies the following Riccati equation:
$$
\nabla_{\grad r} (\mathscr H_r) + \mathscr H_r^2 + R(\cdot,\grad r)\grad r = 0,
$$
(where $R$ is the $(1,3)$-Riemann curvature tensor)
and applying some ODE comparison theory. Plugging these estimates into formula $(*)$ yields estimates for $\mathscr H_f$.
Proofs of most of these facts can be found in Peter Petersen's Riemannian Geometry.
EDIT: Since I wrote this answer, the latest edition of my Riemannian Manifolds book has been published, so let me add that book as another source for the proofs. See Chapter 11, mainly.
A: Another answer might be given in terms of the pushforward of the logarithmic map. In fact, given a curve $t\mapsto\gamma(t)$ on a manifold $\mathbb{M}$, we have that the covariant derivative of $\log_\gamma x$, with $x\in\mathbb{M}$ fixed, reads
\begin{equation}
\nabla_t^\gamma \log_\gamma x = \tfrac{\mathrm{d}}{\mathrm{d}t}\log_\gamma x + \Gamma^{\mathrm{II}}_\gamma(\log_\gamma x,\dot{\gamma}),
\end{equation}
where $\Gamma^{\mathrm{II}}$ denotes the Christoffel symbol of the second kind. Now notice that $\log_{\bullet} x: \mathbb{M}\to T\mathbb{M}$ hence the derivative $\tfrac{\mathrm{d}}{\mathrm{d}t}\log_\gamma x$ may be written in terms of the pushforward endomorphism $\mathrm{d}_y(\log_y x): T_y\mathbb{M}\to T_y\mathbb{M}$ as
\begin{equation}
\tfrac{\mathrm{d}}{\mathrm{d}t}\log_\gamma x = (\mathrm{d}_\gamma\log_\gamma x)(\dot{\gamma}), 
\end{equation}
therefore, ultimately
\begin{equation}
\nabla_t^\gamma \log_\gamma x = (\mathrm{d}_\gamma\log_\gamma x)(\dot{\gamma}) + \Gamma^{\mathrm{II}}_\gamma(\log_\gamma x,\dot{\gamma}).
\end{equation}
This is a general expression that may be particularized after knowing details about the manifold and its metric.
