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.
Best Answer
The answer that you want, namely div $U_g$, is not going to expressible in terms of a geometrically invariant quantity (such as, say, the scalar curvature of $g$) because it depends on the underlying coordinates $x$ in which you have presented the metric. For example, if $\det g(x)$ is constant (which is a coordinate dependent thing), then the divergence of $U_g$ (computed in the $xu$-coordinates, which, I assume, is what you mean by the 'Euclidean divergence') will vanish identically (and conversely, as a matter of fact).
I'll try to explain this in the symplectic formulation, since that's the version I find the clearest, but I think that you can make the translation on your own. You start with a Riemannian metric $g = dx\cdot G(x) dx$, where $G$ is a function on $\mathbb{R}^n$ with values in positive definite symmetric matrices. Let's write $G = F^TF$, where $F$ is invertible (but not necessarily symmetric; you can take $F$ to be the positive definite square root of $G$ if you like, but that's not necessary for my argument). The Lagrangian is then $L = \tfrac12 u\cdot u$, where $u = F(x)\ dx$, regarded as an $\mathbb{R}^n$-valued function on the tangent bundle of $\mathbb{R}^n$ (i.e., an orthonormal coframing of the underlying manifold). Applying the Legendre tranform, the symplectic form on the cotangent bundle pulls back to the tangent bundle to be $$ \Omega = dp\wedge dx = d((F(x)\ u)^T) \wedge dx = (du)^T\ F(x)^T \wedge dx + u^T\ d(F(x)^T)\wedge dx $$ The vector field $U_g$ is the $\Omega$-Hamiltonian vector field associated to $L$, i.e., it satisfies $$ \iota(U_g)\ \Omega = -dL = - u\cdot du $$ (where $\iota(X)$ denotes interior product, what I normally call 'lefthook'). By the standard identity, the divergence of $U_g$ with respect to the Liouville volume form, i.e., $\mu = \tfrac1{n!}\Omega^n$, vanishes identically. Now, you can compare the Liouville volume form with the 'Euclidean' volume form in $xu$-coordinates by noting that, by exterior algebra, we have $$ \tfrac1{n!}\Omega^n = \det(F(x))\ \tfrac1{n!}\bigl((du)^T\wedge dx\bigr)^n $$ It follows that the 'Euclidean' divergence of $U_g$ in the $xu$-coordinates is given by the formula $$ -U_g\bigl(\log|\det(F(x))|\bigr) = -\tfrac12\ U_g\bigl(\log|\det(G(x))|\bigr). $$ Since, as you have already computed, $U_g(x) = F^{-1} u$, it follows that the divergence you want is, as a function, linear in the $u$-coordinates. You might write it as $$ -\tfrac12 \nabla \bigl(\log|\det(G(x))|\bigr)\cdot (F^{-1} u). $$