Let $\theta: M \times M \to \mathbb{R}$ the squared distance function $\theta(x,y)=d(x,y)^{2}$ on complete Riemannian manifold $M$. I would like to calcule the gradient of $d^{2}$, where $d^{2}_{y}(x)=d^{2}(x,y)$, for $x,y\in M$ such that $x \notin Cut(y)\cup \{y\}$. The idea is to consider $\alpha(t,s)$ a variation through geodesics of a minimizing geodesic $\gamma(t)$ with $\gamma(0)=y$ and $\gamma(l)=x$, $l=d(x,y)$ and if $L(s)$ denotes the length of the geodesic $t \to \alpha(t,s)$ , then by first variation formula :
$$\frac{d L(s)}{ds}\Bigr|_{s=0}= \langle V,T\rangle^{l}_{0} - \int\limits^{l}_{0} \langle V, \nabla_{T}T\rangle dt$$
where $T=\frac{\partial \alpha}{\partial t}$ and $V=\frac{\partial \alpha}{\partial s}$. Now , since $\alpha$ os geodesic $\nabla_{T}T=0$, and taking $\alpha$ such that $V$ is the Jacobi Field such that $V(0)=0$ and $V(l)=v$, where $v \in T_{x} M$ , thus :
$$\frac{d L(s)}{ds}\Bigr|_{s=0}=\langle v,T(l)\rangle.$$
Now, I have two questions
$(1)$ Why $T(l)=\frac{-1}{l}exp^{-1}_{x}(y)$?
$(2)$ Why $d(d_{y}(x))(v)=\frac{d L(s)}{ds}\Bigr|_{s=0}$?
Thanks
Just remember the definition of the inverse of exponential map. The geodesic $\alpha(t)$ can be written as $\alpha(t) = \exp_y(tU), ~\|U\| = 1.$ Furthermore, $x = \alpha(l)$. But, what is $\exp_x^{-1}(y)$? It is exactly the velocity of the geodesic connecting $x$ to $y$, i.e, $-\alpha'(l)$ (since $\alpha$ is the geodesic that connects $y$ to $x$). The problem here is that your geodesic finishes at $t = l$ not $t = 1$ so one must normalize it: $\tilde \alpha(t) := \exp_y(\frac{t}{l}V).$ Now one has $\tilde \alpha'(l) = \frac{1}{l}d\exp_y(V) = \frac{1}{l}\alpha'(l).$