Can a Hamiltonian flow be rephrased into a geodesic flow? Let us consider a classical Hamiltonian system composed by $n$ particles, namely, a Hamiltonian function $H:\Lambda\to\mathbb{R}$ where $\Lambda\subset\mathbb{R}^{2n}$ is the phase space. Then, let us write the Hamilton's equations of motion (HEoM)
$$
\dot{q}_{l}=\frac{\partial H}{\partial p_{l}},\qquad \dot{p}_{l}=-\frac{\partial H}{\partial q_{l}}
$$
It is well known that the solutions of the HEoM is a curve on the phase space $x:\mathbb{R}\to\Lambda $ which actually lies on a certain energy level set:
$$
\Sigma(E):=\{x:=(q_{1},\ldots,q_{n},p_{1},\ldots,p_{n})\;|\; H(x)=E\in\mathbb{R}\}.
$$
On the other hand, $\Sigma(E)$ is a hypersurface embedded in $\mathbb{R}^{2n}$, therefore, the Hamiltonian motion is something like a flow on a hypersurface.
Therefore, my question is
Is it possible to deduce the metric tensor $g_{\Sigma(E)}$ on the energy level set $\Sigma(E)$ (having only the Hamiltonian function) such that the HEoM can be rephrased into geodesic equations, i.e.:
$$
\frac{dq^{i}}{ds}=g_{\Sigma(E)}^{ij}p_{j},\qquad \frac{dp_{i}}{ds}=g^{kj}_{\Sigma(E)}\Gamma_{ji}^{l}p_{k}p_{l}
$$
where $\Gamma^{i}_{jk}$ are the Christoffel symbols of the hypersurface?
 A: There is a connection between Hamiltonian mechanics and Riemannian manifolds but it is a little different than your proposal.
Say you have a spatial manifold, parameterized by $x^\mu$, with metric $g_{\mu \nu}$. One can then define the cotangent bundle, where the momenta $p_\mu$ are the cotangent vectors based at points in the base space. Then, the simple Hamiltonian
$$
H = g^{\mu \nu} p_\mu p_\nu
$$
will yield the geodesic equation in $x^\mu$ as its equation of motion. Note that $H$ corresponds to the length of the momentum vector itself $p^2$, and because $\{ H, H \} = 0$, the length of the momentum vector will not change as the particle moves around.
Your idea, that there could be a metric on phase space, would not work for a few reasons. Here is one. Once you specify a metric and a starting point, you also need a velocity to give initial conditions to evolve the geodesic equation-- they are second order differential equations. However, Hamilton's equations are first order, i.e. you only need to specify a starting position in phase space, not a velocity. This is really why we need to double the dimension of the spatial Riemannian manifold, i.e. take its cotangent bundle, in order to get an equivalence between these two different equations.
A: The answer is positive at least for natural Hamiltonian systems (whose Hamiltonian is quadratic in momenta)  and the name for this type of procedure is Eisenhart lift. It is a procedure to describe trajectories of a classical natural Hamiltonian system as geodesics in an enlarged space.
Consider the following Hamiltonian:
$$
H=\frac12 \sum_{i,j=1..n}h^{ij}(q)p_i p_j+e^2 V(q),
$$
Here $h^{ij}$ is a symmetric matrix function depending on generalized coordinates $q$, $V(q)$ is the scalar potential and $e$ is a coupling constant.
Let us introduce additional pair of phase space variables: coordinate $y$ and its conjugate momentum $p_y$ and consider the following Hamiltonian on thus enlarged phase space:
$$
\mathcal{H}=\frac12 \sum_{i,j=1}^{n}h^{ij}(q)p_i p_j+p_y^2 V(q).
$$
It is easy to see that $p_y$ is conserved so for $p_y=e$ the equations of motions for the new Hamiltonian would recover the equations for the original system. At the same time this new Hamiltonian can be interpreted as geodesic Hamiltonian:
$$
 \mathcal{H}=\frac12 \sum_{A,B=1}^{n+1} g^{AB}p_A p_B,
$$
where $p_A=(p_i,p_y)$, and the inverse metric $g^{AB}$ decomposes as $g^{ij}=h^{ij}$, $g^{i, n+1}=0$, $g^{n+1,n+1}=2 V$. Hamiltonian equations for this system are thus equivalent to geodesic equations for the Riemannian $n+1$ dimensional manifold with the length element
$$
ds^2=h_{ij}dq^i dq^j+\frac1{2V}dy^2.
$$
This Hamiltonian is not without problems: if $V$ is zero (changes sign) somewhere, the inverse metric becomes degenerate (changes signature resp.). So, to fix this problem we could add two new coordinates $u,v$ and conjugate momenta $p_u, p_v$ and write the following Hamiltonian:
$$
 \mathcal{H}'=\frac12 \sum_{i,j=1}^{n} h^{ij}p_i p_j+p_v^2 V(q)+p_u p_v.
$$
If we interpret this expression as a geodesic Hamiltonian, then (because of the last term) the $n+2$ dimensional metric has Lorentzian signature with coordinate $u$ being connected to “time” of the original system ($u=-t/e$) and the length element having the form:
$$
 ds^2=h_{ij}dq^idq^j+2du(dv-V du). $$
The equations of motion for the original system are equivalent to null geodesics with $p_v=e$ and $p_u=-H/e$.
Reference:

*

*Cariglia, M., & Alves, F. K. (2015). The Eisenhart lift: a didactical introduction of modern geometrical concepts from Hamiltonian dynamics. European Journal of Physics, 36(2), 025018, doi:10.1088/0143-0807/36/2/025018, arXiv:1503.07802
