Gradient of distance-squared on Riemannian manifold Problem: 
Let $M$ be a Riemannian manifold. Consider the function $f: M \rightarrow \mathbb{R}$ where $f(x)=\text{dist}_M^2(p,x)$, and $p \in M$ is fixed. Show that $\text{grad}(f)=-2\exp^{-1}_x(p)$ as vectors in $T_xM$. (Assuming that $\exp^{-1}$ exists and is smooth etc.)
My attempt at a proof: 
We must show that $\langle -2\exp^{-1}_x(p), \cdot \rangle_x = df(\cdot)$ as 1-forms at $x$.
Let $e_1,\dots,e_n$ be an orthonormal basis of $T_pM$ and introduce normal coordinates at $p$: $x=(x^1,\dots,x^n) \leftrightarrow x=\exp_p(x^i e_i)$. If $q=\exp_p(v)$ then $\text{dist}_M (p,q)=||v||_p$ and so $f(x)=||x^i e_i||^2_p=\sum (x^i)^2$, where the second equality follows from the fact that the metric is the identity at $p$ (in these coordinates).
So $df=\sum 2x^i dx^i$. $(*)$
On the other hand, the geodesic "from x to p" is the same as the geodesic "from p to x", except that the direction is reversed. It should therefore be true that
$p=\exp_x(-S_{p \rightarrow x}x^i e_i )=\exp_x(-x^i S_{p \rightarrow x} e_i)$ where $S_{p \rightarrow x}$ is the parallel transport from $p$ to $x$. Hence $\exp^{-1}_x(p)=-\sum x^i S_{p \rightarrow x} e_i $. Hence
$\langle -2\exp^{-1}_x(p), \cdot \rangle_x=\langle 2x^i S_{p \rightarrow x}e_i, \cdot \rangle = 2x^i \langle S_{p \rightarrow x} e_i, \cdot \rangle = 2x^i dx_i$ which gives the result on comparison with $(*)$. $\square$
My question:
When writing that out, I felt like I was manipulating symbols without really understanding what they mean. For example, when I write $x^i$ I'm not sure whether I mean the coordinate function $x^i$ or its particular value at the point $x$. I was also concerned that when I talk about parallel transport I should really have been talking about the derivative of $\exp$. There's a lemma in Do Carmo which states $\langle (D\exp_p)_v(v),D\exp_p)_v(w)\rangle=\langle v,w \rangle$ for $v \in T_pM$ and $w \in T_v(T_pM) \simeq T_pM$. I can see this might be useful, though I'm not sure precisely how. 
I would appreciate it very much if you checked to see if this proof is correct and "expand it out", particularly the second half) including more details to make clear what is going on.
 A: A straightforward computation one could do is : Let $v\in (TM)_x$. Let $x(t)$ be a curve such that $x(0)=x$ and $x'(0)=v$. Then $df(x)(v)= \frac{d}{dt} d(p,x(t))^2|_{t=0}$. Take a variation of the geodesic $\gamma_0$ joining $p$ to $x$, given by $H(s,t)=\gamma_t(s)$, where $\gamma_t$ is the geodesic joining $p$ to $x(t)$. Notice that $$\frac{d}{dt} d(p,x(t))^2|_{t=0}= 2 \frac{d}{dt} E(\gamma_t)|_{t=0}.$$ Use the first variation formula for energy (see for ex. Gallot Hulin Lafontaine) to find that $$\frac{d}{dt} E(\gamma_t)|_{t=0}=g(v,\gamma_0'(d(p,x))),$$ where $g$ is the metric. But $\gamma_0'(d(p,x))$ is the quantity we wanted to get. It follows now from the definition; $g_x(\nabla(f),v)=df(x)(v)$ for all $v\in (TM)_x$.
A: Here is another proof that uses the Gauss lemma which I believe is what you are referring to in do Carmo.
Let $\gamma_t(s)$ be a geodesic variation through with $\gamma_0(0) = p, \gamma_0(1) = x$ where we note that $\frac{d}{dt}\bigg\vert_{t = 0}\gamma_0(t) = \dot{\gamma}(0) = \log_p(x)$ by definition of $\exp_p$. By reversing the direction of this geodesic we find that $\frac{d}{dt}\bigg\vert_{t = 0}\gamma_0(1-t) = -\dot{\gamma}(1) = \log_x(p)$. This is exactly the statement that you wanted regarding parallel transport, the parallel transport of $\dot{\gamma}(0)$ along $\gamma_0(t)$ gives $\dot{\gamma}(1)$ at $x$ (as the velocity of a geodesic is parallel) which is equal to $-\log_x(p)$.
Take $\frac{d}{ds}\gamma_s(1) = v$. The squared distance function can be written as $f(x) = d(x,p)^2 = \vert \log_p(x) \vert^2$ so that $$df(v) = \frac{d}{ds}\bigg\vert_{s = 0}\vert \log_p(\gamma_s(1))\vert^2  = 2\langle \log_{p*}(v),\log_p(\gamma_0(1))\rangle_p$$
Note that as a function of $s$, $\log_p(\gamma_s(1))$ always lies in $T_pM$ and thus the  inner product in this expression is the inner product of $T_pM$ and so it is as if we are working in $\mathbb{R}^n$ (points in $T_pM$ and tangent vectors can be canonically identified, justifying the last equality above). Under this canonical identification of a vector space and its tangent spaces we have $\log_p(\gamma_0(1)) = \log_{p*}(\dot{\gamma}_0(1))$ which can be seen by finding the coordinate expression for $\gamma_0(t)$ in normal coordinates.
Now we use the Gauss lemma, as $\dot{\gamma}_0(1)$ is a radial the Gauss lemma implies that
$$df(v) =  2\langle \log_{p*}(v),\log_p(\gamma_0(1))\rangle_p = 2\langle \log_{p*}(v),\log_{p*}(\dot{\gamma}_0(1)) \rangle_p = 2\langle v,\dot{\gamma}_0(1) \rangle_x$$
Thus $\text{grad}(f) = 2\dot{\gamma}_0(1) = -2\log_x(p)$ as desired.
Your proof looks to be more or less correct other than the last step which I believe requires some extra justification like the Gauss lemma. The reason being that the normal coordinate vector fields $\partial_i$ are not parallel in general (only the radial ones are!). We thus cannot conclude that the parallel transport of $\partial_{i,p}$ to $x$ is equal to $\partial_{i,x}$.
