  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. Most of this paper is devoted to describing the
method. Details will be provided elsewhere. As examples we revisit the Tremblay
and Winternitz derivation of the Painlev\'e VI potential for a 3rd order
superintegrable flat space system that separates in polar coordinates and 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. The purpose of this project is to develop tools for analysis and
classification of higher order superintegrable systems on any 2D Riemannian
space, not just Euclidean space.
