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.

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