Abstract In the present context, superintegrability is a property of certain probability density functions coming from matrix models, which relates to the average over a distinguished basis of symmetric functions, typically the Jack or Macdonald polynomials. It states that the average can be computed according to a certain combination of those same polynomials, now specialised by specific substitutions when expressed in terms of the power sum basis. For a particular ( q , t ) generalisation of the Gaussian $$\beta $$ β ensemble from random matrix theory, known independently from the consideration of certain integrable gauge theories, we use results developed in a theory of multivariable Al-Salam and Carlitz polynomials based on Macdonald polynomials to prove the superintegrability identity. This then is used to deduce a duality formula for these same averages, which in turn allows for a derivation of a functional equation for the spectral moments.