  Karlin and McGregor's d-variable Hahn polynomials are shown to arise in the
(d+1)-dimensional singular oscillator model as the overlap coefficients between
bases associated to the separation of variables in Cartesian and hyperspherical
coordinates. These polynomials in d discrete variables depend on d+1 real
parameters and are orthogonal with respect to the multidimensional
hypergeometric distribution. The focus is put on the d=2 case for which the
connection with the three-dimensional singular oscillator is used to derive the
main properties of the polynomials: forward/backward shift operators,
orthogonality relation, generating function, recurrence relations,
bispectrality (difference equations) and explicit expression in terms of the
univariate Hahn polynomials. The extension of these results to an arbitrary
number of variables is presented at the end of the paper.
