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.
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.