Derive the Jacobi equation from a specific geodesic variation 
Construction of a specific geodesic variation: Let $\tau = \{\exp_{x_0}(tX):t\in[0,1]\}$ be a geodesic on a complete Riemannian manifold $(M,g)$, starting from $x_0\in M$ with velocity $X\in T_{x_0} M$. Let $V$ be a vector attached to $x_0$ that is not parallel to $X$. Suppose that both $X$ and $V$ have unit length, i.e., $g_{x_0}(X,X) = g_{x_0}(V,V) = 1$. Now we introduce the geodesic starting from $x_0$ with velocity $V$ by
$$y_s := \exp_{x_0}(sV).$$
Denote the parallel transport of $X$ along $\gamma = \{y_s\}$ by
$$X(s) := \Gamma(\gamma)_0^s (X).$$
Then we get a family of geodesics starting from $y_s$ with velocity $X(s)$,
$$\tau^s := \{ \exp_{y_s} (tX(s)):t\in[0,1]\},$$
which forms a variation of the geodesic $\tau$.

So the following vector field along $\tau$ should be the Jacobi field,
$$J(t) := \frac{\partial}{\partial s}\bigg|_{s=0} \exp_{y_s} (tX(s)),$$
which means
\begin{equation}\tag{1}
\frac{D^2}{dt^2} J(t) + R(J(t), \dot \tau(t))\dot \tau(t) =0,
\end{equation}
where $D$ denotes the covariant derivative with respect to the Levi-Civita connection, $R$ the Riemann curvature tensor.

The question is, how to derive the Jacobi equation (1) from this contruction? I tried a lot but failed. This question acturally arises from this paper at its equation (26). Could anyone figure this out? TIA...
 A: This is a rough sketch of the proof, I still have to check for correctness.
We have $\exp_{y_s} (tX(s))$ as a coordinate-field, not a vector field. However, when we do operations like $\frac{D}{ds}$, they correspond to covariant derivatives along the directions of the field given by
$\frac{D}{ds}\exp_{y_s} (tX(s))$.
$$\frac{D^2}{dt} J(t) = \frac{D}{dt} \frac{D}{dt} (\frac{D}{ds}  \exp_{y_s} (tX(s))) = \\=  \frac{D}{dt}  \frac{D}{ds}  \frac{D}{dt}\exp_{y_s} (tX(s)) =  \frac{D}{dt}  \frac{D}{ds} \dot \tau(t) = \\=\frac{D}{ds}  \frac{D}{dt} \dot \tau(t) \; -  R(ds,dt)\dot \tau(t) =\\= -R(J(t), \dot \tau(t))\dot \tau(t)  $$
We can commute the two innermost covariant derivatives because we start off with a coordinate field. (it works roughly for the same reason that $s_{;ab} = s_{;ba}$ for a scalar field). For commuting the two outermost covariant derivatives we apply the Ricci identity.  Also $\frac{D}{dt} \dot \tau(t)$ vanishes because they are geodesics.
The Jacobi equation is a second order differential equation whose solutions are completely determined by the values of $J(p)$ and $\dot J(p)$ for a given parameter $p$ along the geodesic. At $\tau(0) = x_0$  we have $J(0) = V$ and
$$\dot J(0) =  \frac{D}{dt} \frac{D}{ds}  \exp_{y_s} (tX(s))) = \\=   \frac{D}{ds} \frac{D}{dt}  \exp_{y_s} (tX(s)))   = 
\frac{D}{ds} \dot \tau(t)$$
at values $s= 0,t=0$. Given that the initial geodesic family is formed by parallel-transporting $\dot \tau(0)$ along the geodesic spanned by $ds$, we get that $\frac{D}{ds} \dot \tau(t) = 0$ at $s=0,t=0$ . Thus our initial conditions are $J(0) = V$, $\dot J(0) = 0$, which, along with the Jacobi equation, uniquely determine the field along the geodesic.
