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.