  We revise a method by Kalnins, Kress and Miller (2010) for constructing a
canonical form for symmetry operators of arbitrary order for the Schr\"odinger
eigenvalue equation $H\Psi \equiv (\Delta_2 +V)\Psi=E\Psi$ on any 2D Riemannian
manifold, real or complex, that admits a separation of variables in some
orthogonal coordinate system. We apply the method, as an example, to revisit
the Tremblay and Winternitz (2010) derivation of the Painlev\'e VI potential
for a 3rd order superintegrable flat space system that separates in polar
coordinates and, as new results, we give a listing of the possible potentials
on the 2-sphere that separate in spherical coordinates and 2-hyperbolic
(two-sheet) potentials separating in horocyclic coordinates. In particular, we
show that the Painlev\'e VI potential also appears for a 3rd order
superintegrable system on the 2-sphere that separates in spherical coordinates,
as well as a 3rd order superintegrable system on the 2-hyperboloid that
separates in spherical coordinates and one that separates in horocyclic
coordinates. Our aim is to develop tools for analysis and classification of
higher order superintegrable systems on any 2D Riemannian space, not just
Euclidean space.
