Constants of motion for a system of two points moving on spherical surface with a force depending only on their relative distance Consider two points that moves only on a spherical surface of radius $R$. There is only a force between them that has a potential $U(d)$ where $d$ is the distance between the two points.
What is the best choice of coordinates to write the Lagrangian? What are the coordinates to have the greatest number of cyclic coordinates?
And what are the constants of motions of this system? (apart from energy of each point which is a constant of motion clearly since the Lagrangian is time independent) How to find them using the Lagrangian or Hamiltonian?
In spherical coordinate the distance $d$ is
$$d=\|\mathbf{r}-\mathbf{r}^\prime\|=\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}\\
=\sqrt{r^2+r'^2-2rr'\left[\color{red}{\sin(\theta)\sin(\theta')}\color{blue}{\cos(\phi)\cos(\phi')}+\color{red}{\sin(\theta)\sin(\theta')}\color{blue}{\sin(\phi)\sin(\phi')}+\cos(\theta)\cos(\theta')\right]}\\
=\sqrt{r^2+r'^2-2rr'\left[\color{red}{\sin(\theta)\sin(\theta')}\color{blue}{\left(\cos(\phi)\cos(\phi')+\sin(\phi)\sin(\phi')\right)}+\cos(\theta)\cos(\theta')\right]}\\
=\sqrt{r^2+r'^2-2rr'\left[\color{red}{\sin(\theta)\sin(\theta')}\color{blue}{\cos(\phi-\phi')}+\cos(\theta)\cos(\theta')\right]}$$
So if I do a transformation $\Psi=\phi-\phi'$ and $\Psi'=\phi+\phi'$ the potential does not depend on $\Psi'$ and its conjugate momentum should be a constant of motion.
But then I would have only 3 constants of motion.
So is the system a completely integrable hamiltonian system or not?
 A: *

*The extended Lagrangian for OP's system reads
$$ \widetilde{L}~=~ T - U(||{\bf r}_1-{\bf r}_2||)+\sum_{i=1}^2\lambda^i ({\bf r}_i^2 -R^2), \qquad T~:=~\sum_{i=1}^2\frac{m_i}{2} \dot{\bf r}_i^2.\tag{1}$$
The extended model (1) has 2 constraints and 2 Lagrange multipliers $\lambda^i$.

*The extended model (1) is symmetric under time-translations and 3D rotations. By Noether's theorem it therefore has 4 integrals of motion:


*

*The energy $h := T + U$.

*The total angular momentum ${\bf L}:=\sum_{i=1}^2m_i {\bf r}_i\times \dot{\bf r}_i$.


*The corresponding extended Hamiltonian for OP's system reads
$$ \widetilde{H}~=~ T + U(||{\bf r}_1-{\bf r}_2||)-\sum_{i=1}^2\lambda^i ({\bf r}_i^2 -R^2), \qquad T~:=~\sum_{i=1}^2\frac{{\bf p}_i^2}{2m_i} .\tag{2}$$

*We can in principle rewrite the potential term 
$$~U(||{\bf r}_1-{\bf r}_2||)~=~V(\angle({\bf r}_1,{\bf r}_2))\tag{3}$$
as a function of the angle $\angle({\bf r}_1,{\bf r}_2)$ between the positions ${\bf r}_1$ and ${\bf r}_2$. Therefore we can in principle rewrite the systems in terms of coordinates ${\bf q}=(q^1,q^2)$ [e.g. $(\theta, \phi)$] of the 2-sphere $\mathbb{S}^2$ alone without referring to the ambient 3-space $\mathbb{R}^3$. Then the reduced Lagrangian system has $2\times 2=4$ DOF corresponding to 2 point-particles with positions ${\bf q}_1$ and ${\bf q}_2$ on a 2-dimensional manifold. The corresponding reduced Hamiltonian system lives in an 8-dimensional phase space. [We should stress that the reduced Hamiltonian system is constructed from the reduced Lagrangian system, not the extended Hamiltonian system (2).] The reduced  Hamiltonian system is complete/Liouville integrable. By Arnold-Liouville theorem, there exist angle-action variables. The angle variables are (Hamiltonian) cyclic variables.
