I've since found out that where I'm going wrong is that I don't need to find the absolute derivative. I've been told on another physics forum that this problem is framed in terms of Riemann normal coordinates that makes it OK to assume $\frac{D^{2}\xi^{\mu}}{dt{}^{2}}=\frac{d^{2}\xi^{\mu}}{dt{}^{2}}$. Apparently, this is because the distance the cars travel along their separate geodesics on the sphere is a linear function of time $t$.
Anyway, the calculation is as follows.
The distance each car travels north from the equator along its $\phi=$ constant geodesic is $vt$. The angle car-centre of sphere-equator is $vt/r=vt$. Spherical coordinate $\theta=\left(\pi/2\right)-vt$.
The equation of geodesic deviation is $$\frac{d^{2}\xi^{\mu}}{dt{}^{2}}+R_{\phantom{\mu}\beta\alpha\gamma}^{\mu}\xi^{\alpha}\frac{dx^{\beta}}{dt}\frac{dx^{\gamma}}{dt}=0.$$
In effect we want to find $s=\xi^{\phi}$ as a function of $t$.
The non-zero Riemann components for a unit sphere are: $R_{\phantom{\theta}\phi\theta\phi}^{\theta}=\sin^{2}\theta$, $R_{\phantom{\theta}\phi\phi\theta}^{\theta}=-\sin^{2}\theta$,
$R_{\phantom{\theta}\theta\theta\phi}^{\phi}=-1$,
$R_{\phantom{\theta}\theta\phi\theta}^{\phi}=1$.
Let $u^{\sigma}\equiv\frac{dx^{\sigma}}{dt}$.
Then expand out Riemann components for a unit sphere to get:
$$\frac{d^{2}\xi^{\theta}}{dt{}^{2}}=\left(\sin^{2}\theta\right)\left(u^{\phi}u^{\theta}\right)\xi^{\phi}-\left(\sin^{2}\theta\right)\left(u^{\phi}u^{\phi}\right)\xi^{\theta}$$.But as the cars are moving along lines of constant $\phi$,
$$u^{\phi}=\frac{d\phi}{dt}=0$$
and we get$$\frac{d^{2}\xi^{\theta}}{dt{}^{2}}=0$$.
Next set $\mu=\phi$,
expand out the Riemann components to give
$$\frac{d^{2}\xi^{\phi}}{dt{}^{2}}=\xi^{\theta}\left(u^{\theta}u^{\phi}\right)-\xi^{\phi}\left(u^{\theta}u^{\theta}\right).$$
As $\theta=\left(\pi/2\right)-vt$,
$u^{\theta}=\frac{d\theta}{dt}=-v$
and we get$$\frac{d^{2}\xi^{\phi}}{dt{}^{2}}=-v^{2}\xi^{\phi}.$$
Using the wizardry of Wolfram Alpha, the solution to $$\frac{d^{2}y}{dx^{2}}=-ky$$
is$$y=A\sin\left(x\sqrt{k}\right)+B\cos\left(x\sqrt{k}\right).$$
Substituting $v^{2}=k$, $t=x$
and $\xi^{\phi}=y$
gives$$\xi^{\phi}=A\sin\left(vt\right)+B\cos\left(vt\right).$$
To find the constants $A$
and $B$,
we note that when the distance travelled $vt=0$ $\xi^{\phi}=d$,
giving$$\xi^{\phi}=d=A\sin\left(0\right)+B\cos\left(0\right)$$
giving$$B=d,$$
and therefore$$\xi^{\phi}=d\cos\left(vt\right).$$
This is the same answer as calculated in Geodesic devation on a two sphere. In spherical coordinates this is equivalent to $$\xi^{\phi}=d\sin\left(\theta\right).$$
The separation vector is a Jacobi field because it obeys the Jacobi equation.
Here I will derive geodesic deviation from scratch because I find MTW's derivation hard to follow. (Much like everything else in that book.)
Definition 1. Consider a family of timelike geodesics, having the property that in a sufficiently small open region of the Lorentz manifold $(M,g)$ precisely one geodesic passes through every point. Such a collection is called a congruence. The tangent field to this set of curves, parameterized by proper time $s$, is denoted by $u$ and is normalized as $\langle u,u\rangle=-1$.
Let $\gamma(t)$ be some curve transversal to the congruence. This means that the tangent $\dot\gamma$ is never parallel to $u$ in the region under consideration. Imagine that every point on the curve $\gamma(t)$ moves a distance $s$ along the geodesics which passes through that point. Let the resulting point be $H(s,t)$. This defines a map $H:O\longrightarrow M$ for some open $O\subset \mathbb{R}^2$. For each $t$, the curve $s\longmapsto H(s,t)$ is a timelike geodesic with tangent vectors $u(H(s,t))$. We say that $u$ defines a vector field $u\circ H$ along which the map $H$ is tangential, i.e. of the form
$$u\circ H=H_*\circ\frac{\partial}{\partial s}\quad\text{or}\quad u\big(H(s,t)\big)=T_{(s,t)}H\cdot\frac{\partial}{\partial s}$$
In other words, $u$ and $\partial/\partial s$ are $H$-related. We shall use for $u\circ H$ the same letter $u$. Let $v$ be the tangential vector field along $H$ belonging to $\partial/\partial t$
$$v=H_*\circ\frac{\partial }{\partial t}$$
The vectors $v(s,t)$ are tangent to the curves $t\longmapsto H(s,t)$ ($s$ fixed), and represent the separation of points which are moved with the same proper time along the neighboring geodesics of the congruence (beginning at arbitrary starting points). This is shown in the following diagram (taken from Straumann, General Relativity [2013]):
Since $\partial/\partial s$ and $\partial/\partial t$ commute, the tangential fields $u$ and $v$ also commute. To see this, recall that if the vectors of two Lie brackets are $H$-related, then the Lie brackets themselves are $H$-related.
To get the distance between curves, we create a projection operator $id+u\otimes u$ onto the subspace of the tangent space orthogonal to $u$. Thus the relevant infinitesimal separation vector $n$ is
$$
n=v+\langle v,u\rangle u
$$
Note that $\langle n,u\rangle=\langle v,u\rangle-\langle v,u\rangle=0$, so $n$ is indeed perpendicular to $u$. We now show that $n$ is also Lie transported. We have
\begin{align*}
\mathcal{L}_un&=[u,n]=[u,v]+[u,\langle v,u\rangle u]=\big(u\langle v,u\rangle\big)u\\
u\langle v,u\rangle&=\frac{\partial}{\partial s}\langle v,u\rangle
\end{align*}
The normalization condition $\langle u,u\rangle=-1$ implies
$$0=\frac{\partial}{\partial t}\langle u,u\rangle=2\langle \nabla_v u,u\rangle$$
Furthermore, since $\nabla_vu=\nabla_uv$, it follows that
$$\frac{\partial}{\partial s}\langle v,u\rangle=\langle \nabla_uu,v\rangle+\langle u,\nabla_uv\rangle=\langle u,\nabla_vu\rangle=0$$
This shows that indeed
$$\tag{1}
\mathcal{L}_un=0
$$
Next, we consider
$$\nabla^2_uv=\nabla_u\nabla_uv=\nabla_u\nabla_vu=[\nabla_u,\nabla_v]u$$
However, the Riemann tensor for three vectors $X,Y,Z$ is $R(X,Y)Z=[\nabla_X,\nabla_Y]Z-\nabla_{[X,Y]}Z$. Due to $[u,v]=0$, we obtain
$$\nabla^2_uv=R(u,v)u$$
This is called the Jacobi equation for the field $v$.
We now show that $n$ also satisfies this equation. From (1) it follows that
\begin{equation}
\nabla_un=\nabla_uv+\big(u\langle v,u\rangle\big)u+\langle v,u\rangle \nabla_uu=\nabla_uv
\end{equation}
Furthermore,
\begin{equation}
R(u,n)u=R(u,v)u+\langle v,u\rangle R(u,u)u=R(u,v)u
\end{equation}
Putting this all together, we see that $n$ satisfies the Jacobi equation
\begin{equation}
\nabla^2_un=R(u,n)u
\end{equation}
The Jacobi field $n$ is everywhere perpendicular to $u$. In physics we call this the equation for geodesic deviation. For a given $u$ the right hand side defines at each point $p\in M$ a linear map $n\longmapsto R(u,n)u$ of the subspace of $T_pM$ perpendicular to $u$.
In more familiar notation, we have
$$\frac{D^2 n^a}{D\tau^2}=R^a_{\;bcd}u^b n^c u^d$$
So why even bother calling this equation a Jacobi equation? Perhaps there is some general result of Jacobi fields that we wish to apply to GR? This is indeed the case.
Definition 2. Let $\gamma$ be a geodesic. A pair of points $p$, $q\in\gamma$ are said to be conjugate if there exists a Jacobi field $n$ which is not identically zero but vanishes at both $p$ and $q$.
There exists a very powerful theorem regarding conjugate points (suitably generalized here to spacetime, but a similar version holds in positive-definite spaces):
Theorem 1. Let $\gamma$ be a smooth timelike curve connecting two points $p$, $q\in M$. Then the necessary and sufficient condition that $\gamma$ locally maximizes proper between $p$ and $q$ over smooth one parameter variations is that $\gamma$ be a geodesic with no point conjugate to $p$ between $p$ and $q$.
This theorem is the basis for the singularity theorems of GR. Since MTW doesn't discuss these theorems (AFAIK, haven't read the whole thing), they don't need to mention Jacobi fields or conjugate points. Indirectly, this property of Jacobi fields is used in the justification for Big Bang cosmology.
For more in-depth discussions, I refer you to
Wald, General Relativity (1984)
Hawking & Ellis, The Large Scale Structure of Spacetime (1973)
(in that order) and references therein.
Best Answer
The first reason is that your "distance" between geodesics is measured by a parallely propagated direction $\partial/\partial \phi$. If you take a look at the sphere, the difference $\Delta \phi$ does not correspond to the distance between the points on the geodesics. The distance between them would be measured by arc-lengths of great circles. But you are using $\theta=const$ circles which are not the great circles unless $\theta=\pi/2$.
The second reason is you are working on a space of constant curvature. See below.
Say you take a vector $\zeta^\mu$ and propagate it along your geodesic to get $\zeta^\mu(\lambda)$, i.e. you solve $$\frac{D \zeta^\mu}{d \lambda} = 0$$ Thanks to that you now have by a straightforward application of the Leibniz rule $$\frac{D^2(\xi^\mu \zeta_\mu)}{d^2 \lambda} = \frac{D^2\xi^\mu}{d^2 \lambda} \zeta_\mu$$ But $\xi^\mu \zeta_\mu$ is a scalar which is propagated trivially, so you actually also get $D/d\lambda \to d/d\lambda$. Now you can project your equation of geodesic equation into $\zeta^\mu$ to get $$\frac{d^2(\xi^\mu \zeta_\mu)}{d^2 \lambda} + R^\mu_{\;\nu \kappa \lambda} \xi^\kappa u^\nu u^\lambda \zeta_\mu = 0 \; \; \; \;(*)$$
Now to use the fact you are on a space of constant curvature. In such a space, you can express the curvature tensor as $$R_{\mu \nu \kappa \lambda} = K(g_{\mu \kappa} g_{\nu \lambda} - g_{\mu \lambda} g_{\nu \kappa})$$ When you plug this in into the geodesic deviation equation, you get $$\frac{D^2 \xi^\alpha}{d\lambda^2} + u^2 K \xi^\alpha = 0$$ Where $u = dx/d\lambda$ and we choose an orthogonal $\xi$ to it (as is also your case). When you make the projection into the parallely propagated $\zeta_\alpha$, you get $$\frac{d^2 (\xi^\alpha \zeta_\alpha)}{d^2\lambda} + u^2 K \xi^\alpha \zeta_\alpha = 0$$ I.e., if you are investigating only $\Delta \phi = \xi^\alpha \zeta_\alpha$ where $\zeta = \partial/\partial \phi$, you can use even this strictly linear equation.
Note that without using the constant curvature of the space you would end up with equation $(*)$ which doesn't give you much of a clue why you should be able to find $\xi^\alpha \zeta_\alpha$ in the sum with the curvature tensor. So your space is special and your measure of deviation is special - both are necessary ingredients for the assumption.