That's an interesting thought. If we used the geodesic distance instead that could significantly simplify calculations, still not enough to make it exactly solvable but I think it would make perturbation terms very easy to calculate. You could probably exploit the rotational symmetry of the problem and "put" one of the spherical harmonics at the north pole, then the geodesic distance to the other harmonic is simply some multiple of the polar angle. That seems very do-able to arbitray orders in perturbation theory. However the functions you integrate over do need to be eigenstates of the total angular momentum operator so I don't know if that spoils it hm..This is quite an interesting problem as a simplified model of a multielectron atom... Anyone who's more familiar with this, does it make a significant difference in the qualitative features of the energy eigenfunctions, whether I define the electron-electron distance ##r## in the ##e^2 /4\pi\epsilon_0 r## term as the shortest path that stays on the sphere, or as the euclidean distance in the 3D space that the sphere is embedded in?