Higher derivatives than Jacobi fields The first and second derivatives of the distance function (either the full $d:M\times M\to \mathbb{R}$ function or the $d(p,\cdot):M\to \mathbb{R}$ function) as well as the derivative of the exponential map (again both of the full $\exp:TM\to M$ map and of the map $\exp_p:T_pM\to M$) may be computed with the aid of Jacobi fields, i.e, solutions to Jacobi's equation.
I have a scenario where I need second derivatives of the `full' exponential map $\exp:TM\to M$. That is, denoting the pushforward of a differentiable map by a  '$ _*$', I would like to know a thing or two about $\nabla_X\exp_*\mathcal{V}$ (where $\mathcal{V}$ is a section of $TTM$ and $X$ is an appropriate vector field). In particular, I think I will require some comparison techniques analagous to those for Jacobi fields (e.g., Rauch's comparison theorems).
Can anyone point me in the right direction?
 A: The higher derivatives of the exponential map satisfy the corresponding higher derivative of the Jacobi equation (because the first derivative satisfies the Jacobi equation itself), which is just an inhomogeneous Jacobi equation, where the homogeneous part is the original Jacobi equation, and the inhomogeneous term involves lower order covariant derivatives of the Jacobi field and covariant derivatives of the curvature tensor. So you would proceed recursively, bootstrapping pointwise bounds on lower derivatives, as well as pointwise bounds on the curvature tensor and its covariant derivatives, into a pointwise bound of the derivative of the Jacobi field. You'll need to figure out how get pointwise bounds for a solution to an inhomogeneous self-adjoint linear second order ODE. I'm sure this has been done before, probably for exactly the same purpose as here, but I don't know or remember where.
A: You might be interested in the Jacobi flow on $TTM$ whose flow lines project to geodesics, velocity fields of geodesics, and Jacobi fields. You can continue to higher order.


*

*Peter W. Michor: The Jacobi Flow. Rend. Sem. Mat. Univ. Pol. Torino 54, 4 (1996), 365-372
(pdf)
A: Recently in the context of optimization on manifolds, the covariant derivative of the differential of the exponential map $\nabla \text{d } \text{exp}$ was studied in two papers:
1 ''An accelerated first-order method for non-convex optimization on manifolds'' and
2 ''Curvature-Dependant Global Convergence Rates for Optimization on Manifolds of Bounded Geometry''
As @Deane Yang describes, both papers show that this object satisfies a second-order inhomogenous linear ODE, see Proposition B.1 in 1, and Proposition 4.1 in 2.  The following is from 2:

Let $(M, g)$ be a Riemannian manifold, and let $\gamma$ be a geodesic
with initial unit vector $v \in T_p M$.  For two vectors $w_1, w_2$
perpendicular to $v$, the vector field along $\gamma$ $$K(t) = (\nabla
 \text{d } \text{exp}_p)_{tv}(t w_1, t w_2)$$ satisfies the following
ODE along $\gamma$ $$\ddot{K} + R(K, \dot{\gamma})\dot{\gamma} + Y =
 0, \quad K(0) = 0, \dot{K}(0) = 0$$ where $Y$ is the vector field
along $\gamma$ given by $$Y = 2 R(J_1, \dot{\gamma}) \dot{J}_2 +
 2R(J_2, \dot{\gamma})\dot{J}_1 + (\nabla R_{\dot{\gamma}})(J_2,
 \dot{\gamma})J_1 + (\nabla R_{J_2})(J_1, \dot{\gamma})\dot{\gamma}$$
and $J_1, J_2$ are the Jacobi fields along $\gamma$ with initial
conditions $J_i(0) = 0, \dot{J}_i(0) = w_i$.

Using this ODE, both papers derive bounds on $||\nabla \text{d } \text{exp}||$ based on bounds for the sectional curvatures and $\nabla R$.  Paper 1 uses a bootstrapping technique, while paper 2 uses ODE comparison techniques which yield comparatively tighter bounds albeit requiring a little bit more analysis.  For these results, see Theorem B.2 in 1, and Theorems 4.10 and 4.11 in 2.  The following is from 2:

Let $$\text{sn}_{\kappa}(t) = \left\{ \begin{array}{lII}
 \frac{\sin(\sqrt{\kappa}t)}{\sqrt{\kappa}} & \quad \kappa > 0 \\ t & \quad
 \kappa = 0 \\ \frac{\sinh(\sqrt{-\kappa} t)}{\sqrt{-\kappa}} & \quad
 \kappa < 0 \end{array} \right\}$$ and $$\pi_{\kappa} = \left\{
 \begin{array}{lI} \frac{\pi}{\sqrt{\kappa}} & \quad \kappa > 0 \\
 \infty & \quad \kappa \leq 0 \end{array} \right\}.$$
Let $(M, g)$ be a Riemannian manifold with sectional curvatures
bounded below by $\delta$ and above by $\Delta$.  Also assume
$$||(\nabla_u R)(w, u)w + (\nabla_w R)(w, u)u|| \leq 2 \Lambda ||u||^2
 ||w||^2 \quad \forall p \in M \quad u, w \in T_p M.$$
For a geodesic $\gamma : [0, r] \rightarrow M$ with initial unit
vector $v$, $r < \pi_{(\Delta + \delta)/2}$, and any two vectors $w_1,
 w_2 \in T_p M$, we have that $$||(\nabla \text{d } \text{exp}_p)_{r
 v}(w_1, w_2)|| \leq \frac{8}{3r^2} \text{sn}_{\delta}
 \left(\frac{r}{2} \right)^2\left(\Lambda \cdot
 \text{sn}_{\delta}\left(\frac{r}{2}\right)^2 + 2 \max\{|\Delta|,
 |\delta|\}\text{sn}_{\delta}(r) \right) ||w_1|| \cdot ||w_2||.$$

Finally, these bounds are shown to be tight, in a certain sense -- see the discussion after Theorem B.2 in 1, and subsection 4.6.1 in 2.
