On refined Chern-Simons and refined ABJ matrix models

We consider the matrix model of $U(N)$ refined Chern-Simons theory on $S^3$ for the unknot. We derive a $q$-difference operator whose insertion in the matrix integral reproduces an infinite set of Ward identities which we interpret as $q$-Virasoro constraints. The constraints are rewritten as difference equations for the generating function of Wilson loop expectation values which we solve as a recursion for the correlators of the model. The solution is repackaged in the form of superintegrability formulas for Macdonald polynomials. Additionally, we derive an equivalent $q$-difference operator for a similar refinement of ABJ theory and show that the corresponding $q$-Virasoro constraints are equal to those of refined Chern-Simons for a gauge super-group $U(N|M)$. Our equations and solutions are manifestly symmetric under Langlands duality $q\leftrightarrow t^{-1}$ which correctly reproduces 3d Seiberg duality when $q$ is a specific root of unity.
Comment: v2: 30 pages, minor revisions, added comments on relations to quantum mirror curves in Section 2.5, added comments on ABJ integrals in Appendix C, Lett.Math.Phys. version