Let $g$ be a Riemannian metric on the $d$-dimensional flat space $\mathbb R^d$, and consider the usual Lagrangian $$L(x, \dot x) = \tfrac 1 2 g_{ij}(x) \dot x^i \dot x^j.$$ Let $\hat g := \sqrt g$ denote the square root of the metric $g$, implicitly defined by the formula $\hat g_{ai} \hat g_{bj} \delta^{ab} = g_{ij}$, where $\delta^{ab}$ is the identity 2-tensor. I want to introduce phase-space-type coordinates $$u_i = \hat g_{ij} \dot x^j.$$ In the coordinates $(x,u)$, the metric on $u$ is just the Euclidean metric: $\langle u, u \rangle := \delta^{ij} u_i u_j$.
Let $x = x(t)$ denote a geodesic for the metric $g$, and define $u(t) := \hat g_{ij} \dot x^j$. These coordinates are convenient, because along the geodesic, $u(t)$ remains on the sphere of radius $|u(0)|$: $$\tfrac{d}{dt} \langle u, u \rangle = \tfrac{d}{dt} \delta^{ij} u_i u_j = \tfrac{d}{dt} g_{ij} \dot x^i \dot x^j = 0,$$ since geodesics are parametrized by unit speed. Conveniently, this means that $\langle u, \dot u \rangle \equiv 0$.
(You may wonder why don't I just use Hamiltonian phase-space coordinates $(x,p)$. In my research, I consider $g$ as a parameter, ranging over all possible Riemannian metrics on the plane $\mathbb R^2$. Hamiltonian coordinates have the nice property that for a fixed metric, the energy shells $\{ g^{ij}(x) p_i p_j = \mathrm{constant} \}$ are invariant under the geodesic flow. Unfortunately, these energy shells are not independent of the metric parameter $g$. In the coordinates $(x,u)$, on the other hand, the shells $\{ \langle u, u \rangle = \mathrm{constant} \}$ are just spheres in Euclidean space, and do not depend on $g$. In particular, it is important to me that these spherical shells are invariant under rotations in the phase space $\mathbb R^d \times \mathbb R^d$).
I want to calculate the geodesic equation in the coordinates $(x,u)$, particularly for the case that $d=2$. It is easy to see that $\dot x^j = \hat g^{ji} u_i$, where the superscripts denote the inverse of $\hat g$. When I calculate $\dot u$, though, I get a mess: $$\dot u_a = \big( \hat g_{ab,c} \hat g^{cj} \hat g^{bi} - \hat g_{ab} \Gamma_{uv}^b \hat g^{ui} \hat g^{vj} \big) u_i u_j,$$ where $\Gamma_{uv}^b$ are the Christoffel symbols for the metric $g$. I tried simplifying this expression, to no effect. There is plenty of symmetry around (e.g., $\langle u, \dot u \rangle = 0$), and I'm sure that the formula for $\dot u$ takes a much, much simpler form.
Question: Is there a simple expression in these coordinates for the evolution of $u(t)$?
Let me explain why the above expression is inadequate. For the metric $g$, let $U_g$ denote the vector field given by $U_g(x,u) = (u, \dot u)$ (where $\dot u$ is the expression above), so that solutions to the differential equation $(\dot x, \dot u) = U_g(x,u)$ are geodesics for the metric $g$. I need to calculate the (Euclidean) divergence $\operatorname{div} U$. I am pretty sure that in the end, $\operatorname{div} U$ can be expressed in some simple geometric quantities involving the metric (like the Riemannian divergence $\operatorname{div}_g$ of some vector field, scalar curvature $K_g$, etc.). For the messy $\dot u$ above, though, it is impossible for me to see what the true character of $\operatorname{div} U$ is.