Let $(M,g)$ be a Riemannian manifold. Then define
$f: T^*M \times M \to \mathbb{R}$
$f(x,\xi, y) = \langle exp_x^{-1} y, \xi \rangle$
where $exp_{\cdot}\cdot$ is the the exponential map and it's inverse is defined on a neighborhood of the diagonal and $\langle,\rangle$ is the natural pairing between tangent and cotangent vectors.
I'd like to compute the Hessian of $f$ in the x and y variables.
So far I've found that for $\partial^2_{xx} f(x_0,\xi_0,y_0)$ I can put $x$ in normal coordinates at $y_0$ and then $exp_x^{-1} y_0$ is simply $-x$ so
$\partial^2_{xx} f(x_0,\xi_0,y_0) = 0.$
Similarly if $y$ is normal coordinates at $x_0$ I have
$\partial^2_{yy} f(x_0,\xi_0,y_0) = 0.$
But how can I compute the mixed derivatives $\partial^2_{xy}$. I am guessing this will involve Jacobi fields? Thanks.