Consider the following system: \begin{equation} \ddot{x}_1=-\frac{\mu_{\oplus} x_1}{r^3}-\frac{\mu_\oplus R_{\oplus}^2J_2}{r^5}\bigg(\frac{3}{2}x_1-\frac{15}{2}\frac{x_1x_3^2}{r^2}\bigg)\\ \ddot{x}_2=-\frac{\mu_{\oplus} x_2}{r^3}-\frac{\mu_\oplus R_{\oplus}^2J_2}{r^5}\bigg(\frac{3}{2}x_2-\frac{15}{2}\frac{x_2x_3^2}{r^2}\bigg)\\ \ddot{x}_2=-\frac{\mu_{\oplus} x_3}{r^3}-\frac{\mu_\oplus R_{\oplus}^2J_2}{r^5}\bigg(\frac{9}{2}x_3-\frac{15}{2}\frac{x_3^3}{r^2}\bigg). \end{equation} where $J_2=1082.6261\cdot 10^{-6}, \mu_{\oplus}=398600.4405km^3s^{-2},r=\sqrt{x_1^2+x_2^2+x_3^2}, R_{\oplus}=6378.137 km.$
Give me an idea of how I could prove that this system, which describes the case of a satellite whose motion is disturbed by the non-sphericity of the earth, is a Hamiltonian system and how I could determine Hamilton's function in Delaunay coordinates. I understand the practical part well, but from a mathematical point of view I do not have a theoretical basis. Thank you very much !