How to derive the equation for geodesic deviation starting from the geodesic equation? 
A family of geodesics can be parametrised $x^a =  x^a(s, t)$ where $s$ is the distance along a geodesic and the parameter $t$ speciﬁes the geodesic. For each $t$ the geodesic equation is
$$\frac{\partial^2 x^a}{\partial s^2}+\Gamma_{bc}^a\frac{\partial x^b}{\partial s}\frac{\partial x^c}{\partial s}=0\tag{1}$$
Partially diﬀerentiate this equation with respect to $t$ to obtain the equation of geodesic deviation
$$\frac{D^2 w^a}{\partial s^2}=-R_{bcd}^a u^bw^cu^d\tag{2}$$
where $u^a = \partial x^a/\partial s\,$and $w^a=\partial x^a/\partial t$.


So taking the time derivative of $(1)$ gives $$\frac{\partial}{\partial t}\left(\frac{\partial u^a}{\partial s}\right)+\left(\frac{\partial}{\partial t}\Gamma_{bc}^a\right)u^bu^c+\Gamma_{bc}^a\frac{\partial u^b}{\partial t}u^c+\Gamma_{bc}^au^b\frac{\partial u^c}{\partial t}=0\tag{a}$$
But right away I feel totally stuck as I don't understand what to do next in order to reach the desired expression, $(2)$. Could someone please provide hints or tips on how to proceed?

Just for reference, I have typeset the full solution as given by the author below:

Differentiating the geodesic equation with respect to $t$
$$\frac{\partial^2 w^a}{\partial s^2}+\partial_d\Gamma_{bc}^aw^du^bu^c+2\Gamma_{bc}^a\frac{\partial w^b}{\partial s}u^c=0\tag{3}$$
Now
$$\frac{D^2w^a}{\partial s^2}=\frac{\partial}{\partial s}\frac{Dw^a}{\partial s}+\Gamma_{bc}^a\frac{D w^b}{\partial s}u^c$$
$$=\frac{\partial}{\partial s}\left(\frac{\partial w^a}{\partial s}
+\Gamma_{bc}^aw^bu^c\right)+\Gamma_{bc}^a\left(\frac{\partial w^b}{\partial s}+\Gamma_{de}^bw^du^e\right)u^c$$
$$=\frac{\partial^2 w^a}{\partial s^2}+\left(\partial_d\Gamma_{bc}^a\right)u^dw^bu^c+\Gamma_{bc}^aw^b\frac{\partial u^c}{\partial s}+2\Gamma_{bc}^a\frac{\partial w^b}{\partial s}u^c+\Gamma_{bc}^a\Gamma_{de}^bw^du^eu^c$$
$$=\left(\partial_d\Gamma_{bc}^a\right)u^dw^bu^c-\Gamma_{bc}^a\Gamma_{de}^cu^du^ew^b+\Gamma_{bc}^a\Gamma_{de}^bw^du^eu^c+\frac{\partial^2w^a}{\partial s^2}+2\Gamma_{bc}^a\frac{\partial w^b}{\partial s}u^c,$$
using $\partial u^c/\partial s=-\Gamma_{de}^cu^du^e$. Using $(3)$ this can be written as
$$\frac{D^2 w^a}{\partial s^2}=\left(\partial_d\Gamma_{bc}^a\right)u^dw^bu^c-\Gamma_{bc}^a\Gamma_{de}^cu^du^ew^b+\Gamma_{bc}^a\Gamma_{de}^bw^du^eu^c-\partial_d\Gamma_{bc}^aw^du^bu^c.$$
All indices apart from $a$ are dummy indices. In the first two terms on the right-hand side swap $b$ and $c$ indices. In the second two terms swap the $c$ and $d$ indices. This yields
$$\frac{D^2 w^a}{\partial s^2}=\left(\partial_d\Gamma_{bc}^a\right)u^dw^cu^b-\Gamma_{bc}^a\Gamma_{de}^bu^du^ew^c+\Gamma_{bd}^a\Gamma_{ce}^bw^cu^eu^d-\partial_c\Gamma_{bd}^aw^cu^bu^d$$
$$=-\left(\partial_c \Gamma_{bd}^a-\partial_d\Gamma_{bc}^a+\Gamma_{ec}^a\Gamma_{db}^e-\Gamma_{ed}^a\Gamma_{cb}^e\right)u^bw^cu^d=-R_{bcd}^au^bw^cu^d$$

Looking at the solution above, I really don't know how eqn. $(3)$ in the solution was found by differentiating eqn. $(1)$ and moreover, why can I not arrive at the same eqn. via $(\mathrm{a})$?
 A: There are two things you need to use to get from (a) to (3):

*

*Partial derivatives commute:
$$
\frac{\partial}{\partial t} \frac{\partial x^a}{\partial s}
= \frac{\partial}{\partial s} \frac{\partial x^a}{\partial t}
$$

*When a quantity $A$ is a function of $x, y, z$, etc., which in turn all depend on $t$, then
$$
  \frac{\partial}{\partial t} A
  = \frac{\partial A}{\partial x} \frac{\partial x}{\partial t}
    + \frac{\partial A}{\partial y} \frac{\partial y}{\partial t}
    + \frac{\partial A}{\partial z} \frac{\partial z}{\partial t}
    + \mathrm{etc.}
$$
Or, the same idea in slightly different notation: when a quantity $A$ is a function of $x^a$, and $x^a$ depends on $t$, then
$$
  \frac{\partial}{\partial t} A = \partial_a A \frac{\partial x^a}{\partial t}
$$
A: You can get by without coordinates. Let's say that our curvature sign convention is $$R(X,Y)Z = \nabla_X\nabla_YZ-\nabla_Y\nabla_XZ-\nabla_{[X,Y]}Z,$$and let's say that $\gamma\colon I\to M$ is a geodesic. We'll just assume that $\nabla$ is a torsionfree connection, and no metric is required to proceed. Let's write $\nabla_t$ for the covariant derivative operator induced by $\nabla$ on vector fields along $\gamma$. If $\varphi\colon I\times (-\varepsilon,\varepsilon) \to M$ is a variation of $\gamma$, i.e., a smooth map such that $\varphi(t,0) = \gamma(t)$ for all $t\in I$, we may induce partial covariant derivative operators $\nabla_t$ and $\nabla_s$ on vector fields along $\varphi$. Writing $\partial_t$ and $\partial_s$ for partial derivatives, we have the two identities:

*

*$\nabla_t\partial_s\varphi = \nabla_s\partial_t\varphi$ (due to $\nabla$ being torsionfree);

*$\nabla_t\nabla_sX - \nabla_s\nabla_tX = R(\partial_t\varphi,\partial_s\varphi)X$ for any $X \in \Gamma(\varphi^*TM)$.

The geodesic equation is $$\nabla_t\dot{\gamma} = 0.$$
Assuming that every curve $t\mapsto \varphi(t,s)$ is also a geodesic, we have that $$\nabla_t\partial_t\varphi = 0.$$Apply $\nabla_s$ on both sides to get $$\nabla_s\nabla_t\partial_t = 0.$$Now use relation 2. to obtain $$\nabla_t\nabla_s\partial_t\varphi - R(\partial_t\varphi,\partial_s\varphi)\partial_t\varphi = 0.$$Now use relation 1. in the first term above:
$$\nabla_t\nabla_t\partial_s\varphi - R(\partial_t\varphi,\partial_s\varphi)\partial_t\varphi = 0.$$Writing $J_t = \partial_s\varphi(t,0)$, make $s=0$ in the above relation (so that $\partial_t\varphi(t,0) = \dot{\gamma}(t)$) and rearrange terms to obtain the geodesic deviation equation (i.e., the Jacobi equation) $$\nabla_t^2J = R(\dot{\gamma}, J)\dot{\gamma}.$$
If you don't assume that $\nabla$ is torsionfree, you can still do the same, but the final equation will also have a torsion term. People don't care about this since for every connection $\nabla$ there is a unique torsionfree connection with the same geodesics as $\nabla$ (namely, $\nabla - \tau^\nabla/2)$.
A: I agree, the solution given by your source is hard to follow, and, moreover, is non-standard. I will illustrate a more standard approach here.

Consider a spacetime described by a metric $\mathbf g$ with associated connection coefficients $\Gamma$. In this spacetime, imagine a massless particle moving freely through it, having four-position $\boldsymbol x(\tau)$ at the proper time $\tau$ with respect to the particles's wordline. As we know, it obeys the geodesic equation:
$$\ddot x^\alpha +\Gamma^\alpha_{\lambda \kappa}(\boldsymbol x)~ \dot x^\lambda \dot x^\kappa=0\tag{0}$$
Now consider another particle very near to our original particle, separated by the four vector $\boldsymbol d$. (So, the four position of this second particle is given by $\boldsymbol x+\boldsymbol d$.) Because this is another massless particle in our spacetime, it also obeys the geodesic equation. Let's see what this gives us:
$$(\ddot{x}+\ddot{d})^\alpha +\Gamma^\alpha_{\lambda \kappa}(\boldsymbol x+\boldsymbol d)~(\dot{x}+\dot{d})^\lambda(\dot{x}+\dot d )^\kappa=0\tag{1}$$
Before we proceed further, we will expand $\Gamma^\alpha _{\lambda \kappa}(\boldsymbol x+\boldsymbol d)$ about the point $\boldsymbol x$ up to $\mathrm O(|\boldsymbol d|^1)$. Recall Taylor's theorem in several dimensions:
$$F(\boldsymbol x+\boldsymbol \epsilon)=F(\boldsymbol x)+\boldsymbol \epsilon\cdot \nabla F(\boldsymbol x)+\mathrm O(|\boldsymbol \epsilon|^2)$$
Applying this in our case,
$$\Gamma^\alpha _{\lambda \kappa}(\boldsymbol x+\boldsymbol d)=\Gamma^\alpha_{\lambda \kappa}(\boldsymbol x)+d^\mu  \Gamma^\alpha_{\lambda \kappa,\mu}(\boldsymbol x)+\mathrm O(|\boldsymbol d|^2)$$
Where $\Gamma^\alpha _{\lambda \kappa,\mu}:=\partial_\mu\Gamma^\alpha _{\lambda \kappa}$.
Continuing from $(1)$, and discarding terms of $\mathrm O(|\boldsymbol d|^2)$ and $\mathrm O(|\boldsymbol d||\dot{\boldsymbol d}|)$, and dropping the $(\boldsymbol x)$ argument of $\Gamma$ for brevity,
$$ \ddot x ^\alpha+\ddot d^{\alpha}+\big(\Gamma^\alpha_{\lambda \kappa}+d^\mu  \Gamma^\alpha_{\lambda \kappa,\mu}\big)\big(\dot x^\lambda \dot x^{\kappa}+\dot x^\lambda \dot d^\kappa+\dot d^\lambda\dot x^\kappa +\dot d^\lambda \dot d^{\kappa}\big)=0 \\ \ddot x ^\alpha+\ddot d^{\alpha}+\Gamma^\alpha_{\lambda \kappa}\big(\dot x^\lambda \dot x^{\kappa}+\dot x^\lambda \dot d^\kappa+\dot d^\lambda\dot x^\kappa \big)+d^\mu  \Gamma^\alpha_{\lambda \kappa,\mu}\dot x^\lambda \dot x^{\kappa}=0$$
Expanding, and exploiting the symmetry in the lower indices of $\Gamma$,
$$\ddot x ^\alpha+\ddot d^{\alpha}+\Gamma^\alpha_{\lambda \kappa}\dot x^\lambda \dot x^{\kappa}+2\Gamma^\alpha_{\lambda \kappa}\dot x^\lambda \dot d^\kappa+d^\mu  \Gamma^\alpha_{\lambda \kappa,\mu}\dot x^\lambda \dot x^{\kappa}=0$$
Finally, noticing that the first and third terms vanish due to our original geodesic equation, we arrive at our first intermediate result:
$$\ddot d^{\alpha}+2\Gamma^\alpha_{\lambda \kappa}\dot x^\lambda \dot d^\kappa+d^\mu  \Gamma^\alpha_{\lambda \kappa,\mu}\dot x^\lambda \dot x^{\kappa}=0$$
Or
$$\ddot d^{\alpha}+2\Gamma^\alpha_{\lambda \kappa}u^\lambda \dot d^\kappa+d^\mu  \Gamma^\alpha_{\lambda \kappa,\mu}u^\lambda u^{\kappa}=0\tag{2}$$
Where henceforth $\boldsymbol u:=\dot{\boldsymbol x}$.

Let's now define the scalar operator $\mathrm D_{\boldsymbol v}:= v^\alpha \nabla_{\alpha}$, a "directional derivative". It inherits the product rule from the covariant derivative, that is, for any tensors $\mathbf{S,T}$:
$$\mathrm D_{\boldsymbol u}(\mathbf{S}\otimes \mathbf{T})=\mathbf S\otimes \mathrm D_{\boldsymbol u}\mathbf T+\mathbf T\otimes \mathrm D_{\boldsymbol u}\mathbf S$$
You can check this using an example. If we do a direct calculation
$$\mathrm D_{\boldsymbol u}|\boldsymbol v|^2=u^\beta \nabla_{\beta}(v^\alpha v_\alpha) \\ =u^\beta \partial_\beta (v^\alpha v_\alpha) \\ =\frac{\mathrm d}{\mathrm d\tau}(v^\alpha v_\alpha)=\dot{v}^\alpha v_\alpha+v^\alpha\dot{v}_\alpha$$
Where I have used the fact that $\nabla_\mu=\partial_\mu$ for scalars and the fact that $u^\alpha \partial_\alpha=\mathrm d/\mathrm d\tau$ (why?). If you calculate this again using the product rule of the directional derivative, it's certainly more messy,
$$\mathrm D_{\boldsymbol u}(v^\alpha v_\alpha)=v^\alpha \mathrm D_{\boldsymbol u}v_\alpha+v_\alpha \mathrm D_{\boldsymbol u}v^\alpha \\ =v^\alpha \big(u^\beta (\partial_\beta v_\alpha-\Gamma^\gamma_{\beta \alpha}v_\gamma)\big)+v_\alpha\big(u^\beta(\partial_\beta v^\alpha+\Gamma^\alpha_{\beta \gamma}v^\gamma)\big) \\ = v^\alpha\dot{v}_\alpha-\Gamma^\gamma _{\beta \alpha}v^\alpha u^\beta v_\gamma+v_\alpha \dot{v}^\alpha +\Gamma^\alpha_{\beta \gamma}v^\gamma u^\beta v_\alpha \\ =v^\alpha\dot{v}_\alpha+v_\alpha \dot{v}^\alpha-\Gamma^\alpha _{\beta \gamma}v^\gamma u^\beta v_\alpha +\Gamma^\alpha_{\beta \gamma}v^\gamma u^\beta v_\alpha \\ =v^\alpha\dot{v}_\alpha+v_\alpha \dot{v}^\alpha$$
But you get the same result.
Remark: Note the geodesic equation can be written as $\mathrm D_{\boldsymbol u}\boldsymbol u=0$.

Let's calculate $\mathrm D^2_{\boldsymbol u}\boldsymbol d$ and see what happens. In general,
$$\mathrm D_{\boldsymbol u} v^\alpha = \dot{v}^\alpha+\Gamma^\alpha _{\beta \gamma}u^\beta d^\gamma$$
Hence,
$$\mathrm D^2_{\boldsymbol u}d^\alpha=\mathrm D_{\boldsymbol u}(\mathrm D_{\boldsymbol u}d^\alpha) 
\\= \frac{\mathrm d}{\mathrm d\tau}(\mathrm D_{\boldsymbol u}d^\alpha)+\Gamma^\alpha_{\beta \gamma}u^\beta \mathrm D_{\boldsymbol u}d^\gamma
\\=\frac{\mathrm d}{\mathrm d\tau}\big(\dot{d}^\alpha+\Gamma^\alpha_{\beta \gamma}u^\beta d^\gamma\big)+\Gamma^\alpha_{\beta \gamma}u^\beta \big(\dot{d}^\gamma+\Gamma^\gamma _{\delta\epsilon}u^\epsilon d^\delta\big)
$$
We cant just pull $\Gamma$ out of the $\mathrm d/\mathrm d\tau$, because $\Gamma$ depends on $\boldsymbol x$, which depends on $\tau$ in this context. However, we can write
$$\frac{\mathrm d}{\mathrm d\tau}\Gamma^\alpha_{\beta \gamma}=u^\epsilon\partial_\epsilon \Gamma^\alpha_{\beta \gamma}=\Gamma^\alpha_{\beta \gamma,\epsilon}u^\epsilon$$
With this in mind,
$$\mathrm D^2_{\boldsymbol u}d^\alpha= \ddot{d}^\alpha+\frac{\mathrm d}{\mathrm d\tau}\big(\Gamma^\alpha_{\beta \gamma}u^\beta d^\gamma\big)+\Gamma^\alpha_{\beta \gamma}u^\beta \dot{d}^\gamma+\Gamma^\alpha_{\beta \gamma}\Gamma^\gamma_{\delta \epsilon}u^\beta u^\epsilon d^\delta $$
Expanding the derivative and switching some indices around,
$$\mathrm D^2_{\boldsymbol u}d^\alpha=\ddot{d}^\alpha+\Gamma^\alpha_{\beta \delta,\gamma}u^\beta u^\gamma d^\delta+\Gamma^\alpha_{\epsilon \delta}\dot{u}^\epsilon d^\delta+\Gamma^\alpha_{\beta \gamma}u^\beta \dot{d}^\gamma +\Gamma^\alpha_{\beta \epsilon}\Gamma^\epsilon_{\delta \gamma}u^\beta u^\gamma d^\delta$$
Using $(2)$, we can collect the first and fourth term into a single term:
$$\mathrm D^2_{\boldsymbol u}d^\alpha= -\Gamma^\alpha_{\beta\gamma,\delta}u^\beta u^\gamma d^\delta+\Gamma^\alpha_{\beta \delta,\gamma}u^\beta u^\gamma d^\delta +\Gamma^\alpha_{\epsilon\delta }\dot{u}^\epsilon d^\delta +\Gamma^\alpha_{\beta \epsilon}\Gamma^\epsilon_{\delta \gamma}u^\beta u^\gamma d^\delta$$
And, using our original geodesic equation, which can be rewritten as $\dot{u}^\epsilon =-\Gamma^\epsilon_{\beta \gamma}u^\beta u^\gamma$ we can write
$$\mathrm D^2_{\boldsymbol u}d^\alpha = -\Gamma^\alpha_{\beta\gamma,\delta}u^\beta u^\gamma d^\delta+\Gamma^\alpha_{\beta \delta,\gamma}u^\beta u^\gamma d^\delta -\Gamma^\alpha_{\epsilon\delta }\Gamma^\epsilon_{\beta \gamma}u^\beta u^\gamma d^\delta +\Gamma^\alpha_{\beta \epsilon}\Gamma^\epsilon_{\delta \gamma}u^\beta u^\gamma d^\delta$$
Collecting terms, replacing $\epsilon$ with $\mu$,  and switching $\beta$ and $\gamma$ in the fourth term (which is allowed because $\Gamma^\alpha_{\beta \epsilon}\Gamma^\epsilon_{\delta\gamma}u^\beta u^\gamma=\Gamma^\alpha_{\gamma\epsilon}\Gamma^\epsilon_{\delta\beta}u^\beta u^\gamma$)
$$\mathrm D^2_{\boldsymbol u}d^\alpha= \big(\Gamma^\alpha_{\beta \delta,\gamma}-\Gamma^\alpha_{\beta \gamma,\delta}+\Gamma^\mu _{\beta \delta}\Gamma^\alpha_{\mu \gamma}-\Gamma^\mu_{\beta \gamma}\Gamma^{\alpha}_{\mu\delta}\big)u^\beta u^\gamma d^\delta$$
But the thing in the parenthesis is $R^\alpha{}_{\beta\gamma\delta}$! Hence we have derived the geodesic deviation equation:
$$\boxed{\mathrm D^2_{\boldsymbol u}d^\alpha=R^\alpha{}_{\beta \gamma\delta}u^\beta u^\gamma d^\delta}\tag{3}$$

Remark: one can equivalently write
$$\mathrm D_{\boldsymbol u}^2 v^\alpha =\mathrm D_{\boldsymbol u}(\mathrm D_{\boldsymbol u} v^\alpha)=\mathrm D_{\boldsymbol u}(u^\gamma \nabla_\gamma v^\alpha) 
\\ =(\nabla_\gamma v^\alpha)(\mathrm D_{\boldsymbol u}u^\gamma)+u^\gamma \mathrm D_{\boldsymbol u}(\nabla_\gamma v^\alpha) 
\\ =u^\beta u^\gamma\nabla_\beta \nabla_\gamma v^\alpha$$
In other words,
$$\mathrm D_{\boldsymbol u}^2=u^\mu u^\nu\nabla_\mu\nabla_\nu$$
And more generally, it should hold that
$$\mathrm D^n_{\boldsymbol u}=\left(\prod_{i=1}^n u^{\alpha_i}\right)\left(\prod_{i=1}^n \nabla_{\alpha_i}\right)$$
Which is quite appealing.
