UUITP-32/21

On refined Chern-Simons and refined ABJ matrix models
Luca<EMAIL_ADDRESS>and Maxim
<EMAIL_ADDRESS>
Department of Physics and Astronomy, Uppsala University,
Box 516, SE-75120 Uppsala, Sweden

###### Abstract

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.

###### Contents

  1. 1 Introduction
  2. 2 Refined Chern-Simons matrix model
    1. 2.1 The partition function
    2. 2.2 Generating function of correlators
    3. 2.3 q-Virasoro constraints
    4. 2.4 Recursive solution
    5. 2.5 Averages of Macdonald polynomials
    6. 2.6 Symmetries of the constraints
      1. 2.6.1 Langlands duality
      2. 2.6.2 Duality at root of unity
    7. 2.7 Harer-Zagier formulas
  3. 3 Refined ABJ matrix model
    1. 3.1 q-Virasoro constraints
    2. 3.2 Langlands duality
    3. 3.3 Averages of Super-Macdonald polynomials
  4. 4 Conclusions and discussions
  5. A Symmetric functions
  6. B Remarks on torus knots

##  1 Introduction

In the seminal paper [1], Witten gave a gauge field theory construction of the
celebrated Jones polynomials for knots and links in 3 dimensions using the
path integral of Chern-Simons (CS) theory. Since then many generalizations of
this important result have been worked out, both from the mathematical side,
by defining knot (co)homologies that give categorifications of the knot
polynomials, and also from the physics side where people have been able to
compute exactly some BPS observables in certain supersymmetric versions of CS
theory by using localization techniques.

In [2], Aganagic and Shakirov used a deformation of the Verlinde ring of
$\widehat{su}(N)_{k}$ to give a new physical interpretation for certain
refined knot invariants that were previously constructed using homological
methods in [3]. More specifically, using $M$-theory they conjectured a new
type of “refined” gauge theory in 3 dimensions whose Wilson-loop observables
can be identified with the Poincaré polynomials of Khovanov-Rozansky HOMFLY
homology [4, 5]. Their new theory can be seen as a deformation of CS gauge
theory known as refined Chern-Simons (rCS) which depends on two parameters $q$
and $t$ that appear in the definition of the modular $S$ and $T$ matrices. In
the case of usual CS these matrices are written in terms of characters of the
gauge group, which in the case of $U(N)$ are Schur symmetric polynomials. In
the refined theory one substitutes Schur polynomials with their $q,t$
deformation known as Macdonald polynomials which are know to play an important
role in the theory of symmetric functions and in representation theory as
well.

In their paper they also gave a matrix model definition of the partition
function of rCS on $S^{3}$ arguing that (unknot) Wilson-loop observables are
given by averages of Macdonald functions inside of the matrix integral. This
matrix model falls into the category of $q,t$-deformed matrix models for which
it was shown there exists an action of the deformed Virasoro algebra [6] in
the form of Ward identities known as $q$-Virasoro constraints [7].

In this paper we explicitly derive the $q$-Virasoro constraints for the matrix
model of rCS theory on $S^{3}$ by making use of insertions of certain finite
difference operators acting on the integration variables of the matrix
integral. The constraints that we derive can be rewritten as a system of
infinitely many linear relations between the correlators of the model which in
turn can be packaged into some difference equations for the generating
function of the matrix model expectation values.

An explicit analysis of these equations reveals that they can be integrated
completely to give a full solution of the model in terms of normalized
correlation functions, which is equivalent to the statement that the relations
between these correlators form a recursion which is determined uniquely by the
initial data provided by the empty correlator, which is identified with the
partition function itself. Computing the solution up to some finite order we
are able to conjecture a closed formula for the average of Macdonald
characters, and we find that it can be nicely written in terms of the same
characters evaluated at some specific locus. This shows that rCS matrix model
enjoys a property know as superintegrability similarly to what was shown to
happen for other $q,t$-deformed matrix models [8].

Similarly to the case of ordinary unitary gauge groups $U(N)$, one can also
define CS theory for unitary super-groups $U(N|M)$ and in fact there exists a
refinement deformation of this super-group gauge theory as well. This theory
coincides with the refinement of the ABJ theory of [9, 10]. The corresponding
matrix model is therefore the supersymmetric generalization of the rCS one and
we show that it also admits an action of the $q$-Virasoro algebra. Moreover,
we show that the constraints can be written in the same way as those of rCS
once we substitute the rank $N$ with the effective rank $N-\tfrac{\log q}{\log
t}M$. By using this observation, it follows that the constraints in the two
cases have essentially the same solution up to the choice of overall
normalization and rank.

Finally, we show that the refined matrix model representation of the
$q$-Virasoro algebra admits a non-trivial involution which exchanges
$q\leftrightarrow t^{-1}$ known under the name of quantum $q$-geometric
Langlands duality, and it acts non-trivially on the correlators in such a way
that the full (normalized) generating function is left invariant. The physical
interpretation of such a symmetry can be traced back to a 3d version of
Seiberg duality which for pure CS reduces to the well-known level-rank duality
of the representation ring of $\widehat{u}(N)_{k}$

##  2 Refined Chern-Simons matrix model

We start by recalling the definition of the refined CS matrix model as given
in [2, 11]. The refinement is obtained by introducing a specific
$\beta$-deformation of the standard CS matrix model of [12, 13] on the
3-sphere.

###  2.1 The partition function

The theory is described by the following data: the CS level $k\in\mathbb{Z}$,
the gauge group which in our case we take to be $U(N)$ so that the rank and
the dual Coxeter number are both $N$. Additionally one can consider
introducing a Fayet–Iliopoulos (FI) term $\nu$ for the $U(1)$ center of the
gauge group.

The unknot partition function of unrefined CS is given by the Stieltjes–Wigert
matrix integral

$\displaystyle Z_{N}^{\mathrm{CS}}(q)$
$\displaystyle=\frac{(-1)^{N(N-1)/2}}{N!}\int_{\mathbb{R}^{N}}\prod_{1\leq
i<j\leq
N}\left(2\sinh\frac{z_{i}-z_{j}}{2}\right)^{2}\prod_{i=1}^{N}\operatorname{e}^{-\frac{z_{i}^{2}}{2\log
q}}\operatorname{d}z_{i}$ (2.1)
$\displaystyle=\frac{1}{N!}\int_{\mathbb{R}_{+}^{N}}\prod_{i\neq
j}\left(1-x_{i}/x_{j}\right)\prod_{i=1}^{N}\operatorname{e}^{-\frac{\log^{2}x_{i}}{2\log
q}}\frac{\operatorname{d}x_{i}}{x_{i}}$

where in the second line we adopted the change of variables $z_{i}=\log
x_{i}$. The matrix integral depends explicitly on the rank $N$ of the matrices
but also on a complex parameter $q$ which appears in the denominator of the
potential function. In order to identify the matrix model partition function
with the gauge theory partition function one usually assumes that $q$ is of
the form

$q=\operatorname{e}^{\frac{2\pi\operatorname{i}}{k+N}}$ (2.2)

however, we observe that the matrix integral, when convergent, can be regarded
as a well-defined function over the complex plane and not just at roots of
unity.

The refinement of CS theory corresponds to a specific type of deformation of
the matrix integral (2.1), where one substitutes the classical Vandermonde
measure with a $q,t$-analogue of it, i.e. the function

$\Delta_{q,t}(x)=\prod_{i\neq
j}\frac{(x_{i}/x_{j};q)_{\infty}}{(tx_{i}/x_{j};q)_{\infty}}$ (2.3)

where $(x;q)_{\infty}=\prod_{n=0}^{\infty}(1-q^{n}x)$ is the $q$-Pochhammer
symbol. The refined partition function can then be written as

$\displaystyle Z_{N}^{\mathrm{rCS}}(q,t)$
$\displaystyle=\frac{1}{N!}\int_{\mathbb{R}^{N}_{+}}\Delta_{q,t}(x)\prod_{i=1}^{N}\operatorname{e}^{-V(x_{i})}\mathrm{d}x_{i}$
(2.4)

where the potential $V(x)$ is given by the function

$V(x_{i})=\frac{\log^{2}x_{i}}{2\log q}+(-\beta(N-1)-\nu)\log
x_{i}\quad\quad\quad\text{with}\quad\quad\quad\beta=\frac{\log t}{\log q}$
(2.5)

which now also includes an arbitrary FI term proportional to $\log x_{i}$.
When $t=q$ (i.e. $\beta=1$) the function $\Delta_{q,t}(x)$ reduces to the CS
measure $\prod_{i\neq j}(1-x_{i}/x_{j})$ and the partition function in that
limit can be identified with the integral in (2.1) provided we choose the FI
parameter as $\nu=-N$. For later convenience we also introduce the additional
parameters

$r=q^{\nu}\leavevmode\nobreak\ ,\quad\quad\quad p=qt^{-1}\leavevmode\nobreak\
.$ (2.6)

The matrix integral (2.4) depends now on the parameters $q$ and $t$ (and $r$)
which we regard as arbitrary complex numbers. As before, the identification
with the refined gauge theory requires that we choose specific values for
these parameters compatible with the fact that the $S$ and $T$ matrices
provide a unitary representation of $SL(2,\mathbb{Z})$ on the Hilbert space of
a 2-torus. This happens when $q$ and $t$ are given by

$q=\operatorname{e}^{\frac{2\pi\operatorname{i}}{k+\beta
N}}\leavevmode\nobreak\ ,\quad\quad\quad
t=\operatorname{e}^{\frac{2\pi\operatorname{i}\beta}{k+\beta N}}$ (2.7)

which is the refinement of the relation (2.2), however we remark that the
matrix model itself does make sense for arbitrary $q$ and $t$.

######  Remark.

The partition function $Z_{N}^{\mathrm{rCS}}(q,t)$ is a convergent integral
over the positive real line provided $|q|>1$. This can be seen most easily in
the variables $z_{i}=\log x_{i}$ in which case we have a convergent real
Gaussian integral over the contour $\mathcal{C}=\mathbb{R}$. For $|q|\leq 1$
the integrand does not fall-off at infinity on the real line, however one can
make the integral convergent by modifying the integration contour. For
instance, for $|q|=1$ one should rotate the contour slightly off the real axis
while for $|q|<1$ one can take $\mathcal{C}=\operatorname{i}\mathbb{R}$.
Furthermore, one can even substitute the Riemann integral
$\int\operatorname{d}x_{i}$ with the discrete Jackson integral
$\int\operatorname{d}_{q}x_{i}$ without affecting the derivation of the
$q$-Virasoro constraints as discussed in [14], however we do not have
intuition for a physical theory that would give rise to such a “discrete”
matrix model.

The partition function of the refined matrix model computes the vacuum matrix
element of the product $TST$

$Z_{N}^{\mathrm{rCS}}(q,t)=\langle 0|TST|0\rangle$ (2.8)

where $S$ and $T$ are Macdonald deformations of the Kac-Peterson modular
matrices, first introduced in [15, 16]. This matrix element corresponds to the
path integral of rCS theory on $S^{3}$ with no Wilson loops.

As argued in [2, 11], given a representation $R_{\lambda}$ of the gauge group
$U(N)$ one can compute the insertion of the corresponding unknot Wilson loop
operator by inserting the corresponding character in the matrix integral in
(2.4). In the unrefined case, the $U(N)$ characters are given by Schur
polynomials in the integration variables $\\{x_{i}\\}$, while in the refined
case these are substituted by Macdonald symmetric polynomials
$\operatorname{McdP}_{\lambda}(x;q,t)$. The corresponding entry of the modular
$TST$-matrix then computes the Macdonald expectation value

$\langle
R_{\lambda}|TST|0\rangle=\langle\operatorname{McdP}_{\lambda}(x;q,t)\rangle$
(2.9)

Observe that the Wilson loop wraps around the non-contractible cycle of one of
the two solid tori that glue to give the 3-sphere. Moreover, one can interpret
the action of the $T$-matrix as shifting the framing of the unknot by 1 unit.

Similarly, one can compute the rCS path integral of an Hopf link by evaluating
the expectation value of the product of two Macdonald polynomials, however in
this paper we will mostly focus on the unknot partition function.

###  2.2 Generating function of correlators

The correlation functions of the matrix model are the expectation values of
products of power-sums $p_{n}(x)=\sum_{i=1}^{N}x_{i}^{n}$ (see Appendix A),

$c_{\lambda}=\left\langle
p_{\lambda}(x)\right\rangle=\frac{1}{N!}\int_{\mathbb{R}^{N}_{+}}\Delta_{q,t}(x)\sum_{i_{1},\dots,i_{N}=1}^{N}x_{i_{1}}^{\lambda_{1}}\dots
x_{i_{N}}^{\lambda_{N}}\prod_{i=1}^{N}\operatorname{e}^{-V(x_{i})}\mathrm{d}x_{i}$
(2.10)

where $\lambda=[\lambda_{1},\dots,\lambda_{N}]$ is an integer partition of
length at most $N$. It is useful to repackage all the correlators
$c_{\lambda}$ into a generating function. A natural way to do this is to take
the expectation value of the generating function of all power-sums

$\Pi(x,{\bf
t};q,t)=\exp\left(\sum_{n=1}^{\infty}\frac{t_{n}}{n}\frac{1-t^{n}}{1-q^{n}}p_{n}(x)\right)=\sum_{\lambda}z_{\lambda}(q,t)^{-1}p_{\lambda}(x)t_{\lambda}$
(2.11)

with

$t_{\lambda}=\prod_{i=1}^{N}t_{\lambda_{i}}$ (2.12)

and $z_{\lambda}(q,t)$ defined as in (A.13). The function $\Pi(x,{\bf t},q,t)$
is a formal power series in the auxiliary variables ${\bf
t}=\\{t_{1},t_{2},\dots\\}$ which are usually called higher times. If we
regard the higher times $t_{n}$ as the power-sums of a second set of
independent variables, then $\Pi(x,{\bf t},q,t)$ can be interpreted as the
reproducing kernel for the Macdonald inner product $(\cdot,\cdot)_{q,t}$ on
the ring of symmetric functions $\Lambda_{q,t}$.

The generating function of all correlators is defined as

$Z_{N}^{\mathrm{rCS}}({\bf t};q,t)=\left\langle\Pi(x,{\bf
t};q,t)\right\rangle=\sum_{\lambda}z_{\lambda}(q,t)^{-1}c_{\lambda}t_{\lambda}$
(2.13)

Observe that because $\Pi(x,{\bf t};q,t)$ is a reproducing kernel, then the
following identity holds

$p_{n}(x)\Pi(x,{\bf t};q,t)=t_{n}^{\perp}\Pi(x,{\bf t};q,t)$ (2.14)

where $t_{n}^{\perp}$ is the adjoint of multiplication by $t_{n}$, defined by

$t_{n}^{\perp}=n\frac{1-q^{n}}{1-t^{n}}\frac{\partial}{\partial
t_{n}}\leavevmode\nobreak\ ,\quad\quad\quad n\geq 1$ (2.15)

With these definitions we can then write all correlators as

$c_{\lambda}=t_{\lambda}^{\perp}Z_{N}^{\mathrm{rCS}}({\bf t};q,t)|_{{\bf
t}=0}$ (2.16)

Moreover, the operators of multiplication by times $t_{n}$ and their adjoints
$t_{n}^{\perp}$ form an Heisenberg algebra with generators

$a_{n}=\left\\{\begin{array}[]{ll}t^{\perp}_{n}&n\geq 1\\\ N&n=0\\\
t_{-n}&n\leq-1\end{array}\right.$ (2.17)

and commutation relations

$[a_{m},a_{n}]=m\frac{1-q^{|m|}}{1-t^{|m|}}\delta_{m+n,0}$ (2.18)

###  2.3 q-Virasoro constraints

The $q$-Virasoro constraints for the matrix model in (2.4) can be obtained by
introducing an appropriate finite difference operator in the integral of the
generating function. As in the case of the classical Virasoro constraints, the
insertion of the difference operator leads to a vanishing integral even though
the integrand itself is not identically zero. Rewriting the integral as a
linear combination of correlators gives a non-trivial set of relations between
these functions.

In order to define the difference operator appropriate for the rCS matrix
model we first introduce the shift operator

$T_{q,x_{i}}f(x_{1},\dots,x_{i},\dots,x_{N})=f(x_{1},\dots,qx_{i},\dots,x_{N})$
(2.19)

which acts by rescaling the variable $x_{i}$ by the factor $q$. Next we define
the $q$-derivative by

$D_{q,x_{i}}f(x)=\frac{1-T_{q,x_{i}}}{(1-q)x_{i}}f(x)$ (2.20)

and finally we can define our difference operator as

$\mathcal{D}_{m}f(x)=\sum_{i=1}^{N}D_{q,x_{i}}\left[q^{-\frac{m}{2}}x_{i}^{m+1}A_{i}(x;t^{-1})f(x)\right]$
(2.21)

where

$A_{i}(x;t):=\prod_{j\neq
i}\frac{tx_{i}-x_{j}}{x_{i}-x_{j}}\leavevmode\nobreak\ ,$ (2.22)

is an auxiliary function111This is the same function that appears in Chapter
VI of Macdonald’s book [17]. that satisfies

${\Delta_{q,t}(x)^{-1}}{T_{q,x_{i}}\Delta_{q,t}(x)}=t^{-(N-1)}A_{i}(x;t)\left[T_{q,x_{i}}A_{i}(x;t^{-1})\right]^{-1}$
(2.23)

The $q$-Virasoro constraint is obtained by inserting the $q$-difference
operator $\mathcal{D}_{m}$ inside the integral in (2.4). A generalization of
Stokes’ theorem then states that the integral of a total $q$-derivative
vanishes (provided there are no boundary contributions). The resulting
equation can be written as

$\frac{1}{N!}\int\operatorname{d}x^{N}\mathcal{D}_{m}\left\\{\Delta_{q,t}(x)\operatorname{e}^{-\sum_{i}V(x_{i})}\Pi(x,{\bf
t};q,t)\right\\}=\left\langle\mathbb{D}_{m}^{V}\Pi(x,{\bf
t};q,t)\right\rangle=0$ (2.24)

Here we define the operator $\mathbb{D}_{m}^{V}$ by

$\displaystyle\mathbb{D}_{m}^{V}f(x)$
$\displaystyle=\left(\Delta_{q,t}(x)\operatorname{e}^{-\sum_{i}V(x_{i})}\right)^{-1}\mathcal{D}_{m}\left(\Delta_{q,t}(x)\operatorname{e}^{-\sum_{i}V(x_{i})}f(x)\right)$
(2.25)
$\displaystyle=\frac{1}{1-q}\left[\sum_{i=1}^{N}(q^{-\frac{1}{2}}x_{i})^{m}A_{i}(x;t^{-1})-qt^{-(N-1)}\sum_{i=1}^{N}(q^{\frac{1}{2}}x_{i})^{m}\operatorname{e}^{V(x_{i})-V(qx_{i})}A_{i}(x;t)T_{q,x_{i}}\right]f(x)$

This is a discrete family of deformations of the Macdonald operator
$D^{1}_{N}=\sum_{i=1}^{N}A_{i}(x;t)T_{q,x_{i}}$ (also known as the Hamiltonian
of the trigonometric Ruijsenaars-Schneider system), which depends non-
trivially on the choice of potential $V(x)$. As a special case, when $V(x)$ is
$q$-constant then we have

$\mathbb{D}^{V_{q\text{-const}}}_{0}=\frac{t^{-(N-1)}}{1-q}\left[\frac{1-t^{N}}{1-t}-qD^{1}_{N}\right]$
(2.26)

However, the physically interesting case is when the potential is of the form
(2.5) so that

$\operatorname{e}^{V(x_{i})-V(qx_{i})}=rt^{N-1}q^{-\frac{1}{2}}x_{i}^{-1}$
(2.27)

which is equivalent to the statement that the potential function
$\operatorname{e}^{-V(x)}$ can be thought of as a meromorphic section of a
line bundle of degree 1 over the elliptic curve
$\mathbb{C}^{\times}/q^{\mathbb{Z}}$.

If the difference operator $\mathbb{D}^{V}_{m}$ acts on the kernel $\Pi(x,{\bf
t};q,t)$ as multiplication by a symmetric function in the $\\{x_{i}\\}$, then
one can use the identity (2.14) to define a formal adjoint of
$\mathbb{D}^{V}_{m}$. This is indeed possible for $m>0$, in which case we
define the formal adjoint of $\mathbb{D}^{V}_{m}$ as the differential operator
in times given by222For $m\leq 0$, $\mathbb{D}^{V}_{m}\Pi(x,{\bf t};q,t)$ is
not polynomial in the $x_{i}$ but also contains negative powers. For this
reason one cannot take the adjoint of those operators.

$\displaystyle(\mathbb{D}^{V}_{m})^{\perp}=\oint_{|z|=1}\frac{\operatorname{d}z}{2\pi\operatorname{i}z}z^{-m}\frac{1}{(1-q)(1-t^{-1})}\Bigg{\\{}1-t^{-N}\exp\left(\sum_{n=1}^{\infty}\frac{(q^{-\frac{1}{2}}z)^{n}}{n}(1-t^{n})t_{n}^{\perp}\right)$
(2.28)
$\displaystyle+prz\left[1-t^{N}\exp\left(-\sum_{n=1}^{\infty}\frac{(q^{\frac{1}{2}}z)^{-n}}{n}(1-t^{n})t_{n}\right)\exp\left(\sum_{n=1}^{\infty}\frac{(q^{\frac{1}{2}}z)^{n}}{n}(1-t^{-n})t_{n}^{\perp}\right)\right]\Bigg{\\}}$

The $q$-Virasoro constraints then can be rewritten as a set of linear
differential equations for the generating function $Z^{\mathrm{rCS}}_{N}({\bf
t};q,t)$, namely

$(\mathbb{D}_{m}^{V})^{\perp}Z_{N}^{\mathrm{rCS}}({\bf t};q,t)=0$ (2.29)

for all $m>0$.

Observe that the factor of $z$ appearing in the beginning of the second line
of (2.28) is due to the variation of the CS potential (2.27) and it induces a
shift in degree in the power series expansion in $z$. More explicitly, after
taking the residue in $z$, the operator $\mathbb{D}^{V}_{m}$ becomes the sum
of two pieces: the first is a differential operator of degree $m$ while the
second is of degree $m-1$ in the higher times. This shift in degree is
precisely what will allow us in the next section to establish a recursion on
the correlators $c_{\lambda}$.

###  2.4 Recursive solution

In this section we set out to solve the constraints in (2.29). The first step
is to evaluate explicitly the contour integral in (2.28),

$\displaystyle(\mathbb{D}_{m}^{V})^{\perp}=-t^{-N}q^{-\frac{m}{2}}h_{m}$
$\displaystyle\left(\left\\{p_{n}=(1-t^{n})t_{n}^{\perp}\right\\}\right)+pr\left(-pq^{-\frac{1}{2}}\right)^{m-1}\Bigg{[}\delta_{m,1}+$
(2.30)
$\displaystyle-\sum_{\ell=0}^{\infty}t^{N-\ell}e_{\ell}\left(\left\\{p_{n}=(1-t^{n})t_{n}\right\\}\right)e_{\ell+m-1}\left(\left\\{p_{n}=(1-t^{n})t_{n}^{\perp}\right\\}\right)\Bigg{]}$

where $h_{m}$ and $e_{m}$ are the complete homogeneous and elementary
symmetric functions of degree $m$, respectively, which we think of as
polynomials in the power-sums $p_{n}$ as in (A.9).

Assuming the series expansion of $Z_{N}^{\mathrm{rCS}}({\bf t};q,t)$ as in
(2.13), we can let all the differential operators $t^{\perp}_{n}$ act on the
monomials in times and the constraint becomes a linear system for the
coefficients $c_{\lambda}$. In particular this system is triangular w.r.t. a
certain ordering on the partitions and the solution is unique up to overall
normalization specified by the initial data $c_{\emptyset}:=\langle
1\rangle=Z_{N}^{\mathrm{rCS}}(0;q,t)$, i.e. the empty correlator. The proof
goes along the lines of the one described in [8].

Solving the recursion explicitly in degree 1 and 2 we get the following values
for the correlators

$\displaystyle c_{[2]}$
$\displaystyle=\frac{\left(t^{N-1}(qt+q+1)-1\right)}{1+t}\frac{1-t^{N}}{1-t}(rq^{\frac{3}{2}}t^{N-1})^{2}\,c_{\emptyset}$
(2.31) $\displaystyle c_{[1,1]}$
$\displaystyle=\frac{\left(1-t^{N-1}(q(t-1)+1)\right)}{1-t}\frac{1-t^{N}}{1-t}(rq^{\frac{3}{2}}t^{N-1})^{2}\,c_{\emptyset}$
$\displaystyle c_{[1]}$
$\displaystyle=\frac{1-t^{N}}{1-t}(rq^{\frac{3}{2}}t^{N-1})\,c_{\emptyset}$

Similarly one can plug (2.30) into a computer algebra system and compute all
correlators $c_{\lambda}$ up to arbitrarily high degree of the partition
$\lambda$. While this solution is fully explicit and straightforward to
compute, it is not very illuminating. In section 2.5 we discuss a more
systematic way to express the solution to the constraints.

###  2.5 Averages of Macdonald polynomials

By using the explicit solution derived in 2.4 we can now ask what is the
expectation value of a Macdonald polynomial in representation $R_{\lambda}$.
Using the fact that Macdonald symmetric functions admit an expansion into
power-sums we can easily rewrite any Macdonald expectation value as a linear
combination of the correlators $c_{\lambda}$ that we already computed.
Computer experiments then suggest that the following formula should hold

$\frac{\langle\operatorname{McdP}_{\lambda}(x)\rangle}{\langle
1\rangle}=\frac{\operatorname{McdP}_{\lambda}\left(\left\\{p_{n}=\frac{-(-pq^{\frac{1}{2}}rt^{N})^{n}}{1-t^{n}}\right\\}\right)}{\operatorname{McdP}_{\lambda}\left(\left\\{p_{n}=\frac{1}{1-t^{n}}\right\\}\right)}{\operatorname{McdP}_{\lambda}\left(\left\\{p_{n}=\frac{1-t^{nN}}{1-t^{n}}\right\\}\right)}$
(2.32)

which we checked for all partitions up to $|\lambda|\leq 6$. This is a special
case of a property of certain matrix models called superintegrability first
observed in [18, 19, 20, 21]. Superintegrability for $q,t$-deformed matrix
models was later conjectured in [22] and in [8] its was explicitly checked (up
to some finite order) using the solution of the corresponding $q$-Virasoro
constraints.

The last factor on the right of (2.32) is the well-known Macdonald dimension
of the representation $R_{\lambda}$,

$\dim_{q,t}(R_{\lambda})=\operatorname{McdP}_{\lambda}\left(\left\\{p_{n}=\frac{1-t^{nN}}{1-t^{n}}\right\\}\right)=t^{\frac{1}{2}|\lambda|(N-1)}\operatorname{McdP}_{\lambda}\left(x_{i}=t^{\rho_{i}}\right)$
(2.33)

where

$\rho_{i}=\frac{1}{2}(N-2i+1)$ (2.34)

is the Weyl vector of $U(N)$ and for $\beta\in\mathbb{N}$ we can write

$\operatorname{McdP}_{\lambda}\left(t^{\rho}\right)=t^{-\lambda\cdot\rho}\prod_{\alpha>0}\frac{(t^{\alpha\cdot\rho+1};q)_{\alpha\cdot\lambda}}{(t^{\alpha\cdot\rho};q)_{\alpha\cdot\lambda}}=\prod_{\alpha>0}\prod_{m=0}^{{\beta}-1}\frac{\sin\left(\frac{\pi(m+\alpha\cdot(\beta\rho+\lambda))}{k+\beta
N}\right)}{\sin\left(\frac{\pi(m+\alpha\cdot(\beta\rho))}{k+\beta N}\right)}$
(2.35)

where $\alpha>0$ are the positive roots of $U(N)$.

The fractional term on the right of (2.32) is instead a power of $q$ which for
$q$ a root of unity is just a complex phase. More precisely, a direct
computation shows that (2.32) becomes

$\frac{\langle\operatorname{McdP}_{\lambda}(x)\rangle}{\langle
1\rangle}=\operatorname{McdP}_{\lambda}\left(t^{\rho}\right)\left[q^{\frac{1}{2}\lambda\cdot\lambda}t^{\lambda\cdot\rho}(t^{N}rp)^{|\lambda|}\right]$
(2.36)

Fixing the FI parameter as $r=t^{-N}p^{-1}$ we precisely obtain that the
overall phase corresponds to a framing factor of 1 unit

$\frac{T_{\lambda}}{T_{0}}=q^{\frac{1}{2}\lambda\cdot\lambda}t^{\lambda\cdot\rho}$
(2.37)

where $T_{\lambda}$ is the eigenvalue of the modular $T$-matrix corresponding
to the representation $R_{\lambda}$.

Fixing instead $r=t^{-N}q^{-1}$ we get the identity

$\displaystyle\operatorname{McdP}_{\lambda}\left(t^{\rho}\right)\left[q^{\frac{1}{2}\lambda\cdot\lambda}t^{\lambda\cdot\rho-|\lambda|}\right]$
$\displaystyle=\operatorname{McdP}_{\lambda}\left(t^{{\boldsymbol{\delta}}}\right)\left[q^{\frac{1}{2}\lambda\cdot\lambda}t^{-\frac{|\lambda|}{2}-\frac{1}{2}\lambda^{\prime}\cdot\lambda^{\prime}}\right]$
(2.38)

where the r.h.s. now matches with the result of [2] (${\boldsymbol{\delta}}$
is the staircase partition of (A.4)).

Because Macdonald polynomial are a complete basis of the ring of symmetric
functions in the $\\{x_{i}\\}$, we can use formula (2.32) to compute
expectation values of all possible observables of the model. In particular, we
can give a full solution for the generating function by expanding in Macdonald
polynomials (as opposed to power-sum polynomials) as

$Z_{N}^{\mathrm{rCS}}({\bf
t};q,t)=\sum_{\lambda}\left\langle\operatorname{McdP}_{\lambda}(x)\right\rangle\operatorname{McdQ}_{\lambda}({\bf
t})$ (2.39)

This is usually referred to as character expansion of the generating function.

###  2.6 Symmetries of the constraints

We now consider the operators (2.28) and study their symmetries which will
imply symmetries of the solution.

####  2.6.1 Langlands duality

It is well known that $q,t$-deformed models enjoy a non-perturbative symmetry
associated to the exchange of $q$ and $t^{-1}$ which is known as “quantum
$q$-geometric” Langlands duality [23, 24, 25]. To show that our $q$-Virasoro
constraints satisfy this symmetry we define the involution
${{}^{\mathrm{L}}\omega}$ as

$\begin{array}[]{ccccc}{{}^{\mathrm{L}}\omega}(\beta)=\tfrac{1}{\beta}&&{{}^{\mathrm{L}}\omega}(N)=-\beta
N&&{{}^{\mathrm{L}}\omega}(\nu)=-\frac{1}{\beta}\nu\end{array}$ (2.40)

and one can check that $({{}^{\mathrm{L}}\omega})^{2}=\operatorname{Id}$. We
can extend the action of the involution to the full generating function by
imposing that the times transform as well, with the rule

${{}^{\mathrm{L}}\omega}(t_{n})=-\frac{1-t^{n}}{1-q^{n}}p^{\frac{n}{2}}t_{n}$
(2.41)

so that the $q$-Virasoro operators become invariant under Langlands duality,

${{}^{\mathrm{L}}\omega}\left((\mathbb{D}^{V}_{m})^{\perp}\right)=(\mathbb{D}^{V}_{m})^{\perp}$
(2.42)

This implies that the normalized solution to the constraints must also be
invariant under the duality333Observe that this is not true for the partition
function, in fact  ${{}^{\mathrm{L}}\omega}(Z^{\mathrm{rCS}}_{N}(0;q,t))\neq
Z^{\mathrm{rCS}}_{N}(0;q,t)$ . More precisely, we can write

$\frac{Z_{N}^{\mathrm{rCS}}(\\{t_{n}\\};q,t)}{Z_{N}^{\mathrm{rCS}}(\\{0\\};q,t)}=\frac{Z^{\mathrm{rCS}}_{-\beta
N}(\\{-\frac{1-t^{n}}{1-q^{n}}p^{\frac{n}{2}}t_{n}\\};t^{-1},q^{-1})}{Z^{\mathrm{rCS}}_{-\beta
N}(\\{0\\};t^{-1},q^{-1})}$ (2.43)

Observe that for arbitrary complex values of $\beta$, the dual theory has rank
$-\beta N$ which generically is not a positive integer. Clearly the matrix
model of this theory is not well-defined since the rank should correspond to
the number of eigenvalues $x_{i}$. Nevertheless, the $q$-Virasoro constraints
are well-defined for any complex value of the rank, hence the generating
function in the r.h.s. of (2.43) does make sense as a formal power series
solution to the constraints. In other words, one can compute the correlators
for integer values of the rank and then analytically continue their expression
over the whole complex plane.

######  Remark.

As a special case of the identity (2.43), we can choose the $\beta$-parameter
in such a way that the rank of the model on the right becomes 1. Namely, we
can fix the value of $\beta$ as

$\beta=-\frac{1}{N}$ (2.44)

so that we obtain the “universal” relation

$\frac{Z_{N}^{\mathrm{rCS}}(\\{t_{n}\\};t^{-N},t)}{Z_{N}^{\mathrm{rCS}}(\\{0\\};t^{-N},t)}=\frac{Z_{1}^{\mathrm{rCS}}(\\{\frac{1-t^{n}}{1-t^{nN}}t^{\frac{n}{2}(N-1)}t_{n}\\};t^{-1},t^{N})}{Z_{1}^{\mathrm{rCS}}(\\{0\\};t^{-1},t^{N})}$
(2.45)

Moreover, in rank $N=1$ the solution is particularly simple,

$Z_{1}^{\mathrm{rCS}}({\bf
t};q,t)=\left\langle\exp\left(\sum_{n=1}^{\infty}\frac{x^{n}}{n}\frac{1-t^{n}}{1-q^{n}}t_{n}\right)\right\rangle=\sum_{m=0}^{\infty}\langle
x^{m}\rangle
h_{m}\left(\left\\{p_{n}=\frac{1-t^{n}}{1-q^{n}}t_{n}\right\\}\right)$ (2.46)

where we can evaluate the momentum $\langle x^{m}\rangle$ explicitly by the
integral

$\langle
x^{m}\rangle=\int_{0}^{\infty}\mathrm{d}x\,x^{m+\nu}\exp\left(-\frac{\log^{2}x}{2\log
q}\right)=q^{\frac{1}{2}(m+\nu+1)^{2}}\sqrt{2\pi\log q}$ (2.47)

####  2.6.2 Duality at root of unity

Another interpretation of the symmetry of the operators
$(\mathbb{D}^{V}_{m})^{\perp}$ comes from restricting to the physical region
of parameters where $q,t$ are roots of unity as defined in (2.7). It is known
that the $q$-Virasoro algebra at root of unity degenerates to a simpler
algebra related to parafermions [26, 27, 28, 29]. We do not pursue this avenue
here however we remark that in this region of parameters the Langlands duality
reduces to 3d Seiberg duality for pure CS theories, which also coincides with
level-rank duality (see [30] for a discussion of dualities in unrefined CS).
The identification holds provided we assume the identity

$q^{k}t^{N}=1$ (2.48)

which follows from (2.7). For this choice of parameters we have

$q^{k}=q^{-\beta N}$ (2.49)

so that formally we can modify the involution to act as

${{}^{\mathrm{L}}\omega}(N)=k$ (2.50)

therefore, Langlands duality maps refined CS theories as

$\begin{array}[]{ccc}U(N)_{k}&\leftrightarrow&U(k)_{N}\\\ &&\\\
Z^{\mathrm{rCS}}_{N}(q,t)&\simeq&Z^{\mathrm{rCS}}_{k}(t^{-1},q^{-1})\end{array}$
(2.51)

From (2.48) it also follows that rCS theory has an additional discrete
$\mathbb{Z}$-symmetry corresponding to shifts of the rank $N$ of the form

$N\mapsto N+\delta N$ (2.52)

with

$\delta N\in\tfrac{1}{\beta}\left(k+\beta N\right)\mathbb{Z}$ (2.53)

It would be interesting to find a physical motivation that might explain this
discrete symmetry, perhaps via some string/M-theory construction.

###  2.7 Harer-Zagier formulas

Matrix models that possess the superintegrability property described in the
previous section often usually have simple formulas that express the
expectation value of single trace operators. In our language, these are the
correlators whose partition is of length one, i.e. those of the form
$c_{[n]}=\langle p_{n}\rangle$. Knowledge of a closed formula for such
correlators can be very important because they are the coefficients of the
series expansion of the 1-point resolvent of the model. It is believed that in
models that exhibit superintegrability the (Laplace transforms of) single-
trace correlators $c_{[n]}$ are rational functions with zeroes and poles at
powers of $q$ [31, 32], and these formulas are usually called Harer-Zagier
(HZ) formulas [33].

While giving an analytic derivation of such HZ formulas can be quite
complicated, guessing the actual form of the result can be done more easily.
In this section we conjecture a formula of HZ-type for the CS matrix model
when the refinement parameter $\beta$ is set to 1 (i.e. $t=q$). First we
conjecture the following closed formula for 1-point correlators (see also
Proposition 1.1 of [34])

$c_{[n]}=-(-q^{N+\frac{1}{2}}r)^{n}\frac{1-q^{N}}{1-q^{n}}\,\mathcal{P}^{(0,1)}_{n-1}(q^{N-1};q)\,c_{\emptyset}$
(2.54)

where $\mathcal{P}_{n}^{(\alpha,\beta)}(z;q)$ is a little $q$-Jacobi
polynomial defined as

$\mathcal{P}^{(\alpha,\beta)}_{n}(z;q)={}_{2}\phi_{1}(q^{-n},q^{\alpha+\beta+n+1};q^{\alpha+1};q,qz)=\sum_{m=0}^{n}\frac{(q^{-n};q)_{m}(q^{\alpha+\beta+n+1};q)_{m}}{(q^{\alpha+1};q)_{m}(q;q)_{m}}(qz)^{m}$
(2.55)

It is then useful to consider the discrete Laplace transform of this formula
w.r.t. the rank $N$. The transform takes the form of a generating
function444One should regard this formula as a special case of the HZ formula
for torus knots in Section 6 of [32].

$\sum_{N=0}^{\infty}\lambda^{N}c_{[n]}=r^{n}q^{\frac{n(n+2)}{2}}\frac{\lambda}{1-\lambda}\frac{(\lambda;q)_{n}^{2}}{(\lambda;q)_{2n+1}}\,c_{\emptyset}$
(2.56)

where now we can observe that the r.h.s. is indeed a rational function of
$\lambda$ with zeroes and poles at integer powers of $q$. It would be
interesting to derive a similar formula for the refined correlators at $t\neq
q$.

##  3 Refined ABJ matrix model

The ABJ matrix model of [10, 35, 36] admits a similar refinement deformation
as that of the CS matrix model. The partition function of ABJ theory can be
expressed as a matrix integral over two sets of independent variables
$\\{x_{i}\\}_{i=1}^{N}$ and $\\{y_{a}\\}_{a=1}^{M}$ which can be interpreted
as the eigenvalues of an Hermitean super-matrix in the Lie super-algebra of
$U(N|M)$. The usual Vandermonde determinant in the measure is first
substituted by the Cauchy determinant and then $q,t$-deformed to the function
(see [37, 38, 39])

$\Delta_{q,t}(x,y)=\frac{\Delta_{q,t}(x)\Delta_{t^{-1},q^{-1}}(y)}{\prod_{i=1}^{N}\prod_{a=1}^{M}(1-\sqrt{t/q}\leavevmode\nobreak\
x_{i}/y_{a})(1-\sqrt{t/q}\leavevmode\nobreak\ y_{a}/x_{i})}$ (3.1)

The potential of the model is also modified in order to account for the second
set of $\\{y_{a}\\}$ variables and we write it as

$\displaystyle V(x,y)=$ $\displaystyle\sum_{i=1}^{N}\frac{\log^{2}x_{i}}{2\log
q}+\big{(}M-\beta(N-1)-\nu\big{)}\sum_{i=1}^{N}\log x_{i}$ (3.2)
$\displaystyle+\sum_{a=1}^{M}\frac{\log^{2}y_{a}}{2\log
t^{-1}}+\big{(}N-\tfrac{1}{\beta}(M-1)+\tfrac{1}{\beta}\nu\big{)}\sum_{a=1}^{M}\log
y_{a}$

so that the first set of eigenvalues $x_{i}$ has coupling $\log q$ while the
eigenvalues $y_{a}$ have coupling $\log t^{-1}$ (the FI parameters are also
related in such a way as to preserve symmetry under Langlands duality).
Observe that this is the refinement deformation of the usual statement that,
if we regard ABJ theory as a quiver gauge theory with two nodes $U(N)$ and
$U(M)$, then the two CS levels must be opposite to each other555In the case of
the ABJ theory, supersymmetry prevents the CS level to receive quantum
corrections from the dual Coxeter number. Therefore one has
$q=\operatorname{e}^{{2\pi\operatorname{i}}/{k}}$ instead of (2.7). A similar
conclusion was reached in [40] in the case of super-group CS theory..

Expectation values of gauge invariant operators are given by averages of
symmetric functions in the eigenvalues $\\{x_{i}\\}$ and $\\{y_{a}\\}$. In the
following we restrict ourselves to those symmetric functions which can be
written as polynomials in the deformed power-sums of [38] (see also [41]),
defined as

$p_{n}(x,y):=p_{n}(x)-\frac{1-q^{n}}{1-t^{n}}p^{-\frac{n}{2}}p_{n}(y)$ (3.3)

The reproducing kernel which gives the generating function of such deformed
power-sums is the function

$\Pi(x,y,{\bf
t};q,t)=\exp\left(\sum_{n=1}^{\infty}\frac{t_{n}}{n}\frac{1-t^{n}}{1-q^{n}}p_{n}(x,y)\right)$
(3.4)

which can be obtained from (2.11) by substituting $p_{n}(x)$ with
$p_{n}(x,y)$. It follows that

$p_{n}(x,y)\Pi(x,y,{\bf t};q,t)=t^{\perp}_{n}\Pi(x,y,{\bf t};q,t)$ (3.5)

Finally, we can write down the generating function of the rABJ matrix model as

$Z_{N,M}^{\mathrm{rABJ}}({\bf
t};q,t)=\frac{1}{N!M!}\int_{\mathbb{R}^{N+M}_{+}}\Delta_{q,t}(x,y)\Pi(x,y,{\bf
t};q,t)\operatorname{e}^{-V(x,y)}\prod_{i=1}^{N}\mathrm{d}x_{i}\prod_{a=1}^{M}\mathrm{d}y_{a}\,$
(3.6)

As already mentioned, in the unrefined limit $t=q$ the measure
$\Delta_{q,t}(x,y)$ reduces to a Cauchy determinant and one ends up with the
usual matrix model for the ABJ theory which is the super-group version of
(2.1) (see [36] for details). Moreover, it is evident from the definitions
that the rABJ matrix model contains rCS as the special case $M=0$, i.e.

$Z^{\mathrm{rCS}}_{N}({\bf t};q,t)=Z^{\mathrm{rABJ}}_{N,0}({\bf t};q,t)$ (3.7)

######  Remark.

Convergence of the matrix integral for the rABJ partition function requires a
careful choice of contour of integration depending on the values of the
parameters $q$ and $t$. In the following we assume that such a choice of
contour always exists, however we do not write down the definition of the
contour explicitly as it is known that $q$-Virasoro constraints are
independent of such choice. Moreover, as for the case of rCS, one could even
define a discrete version of the rABJ model using the substitutions
$\int\operatorname{d}x_{i}\mapsto\int\operatorname{d}_{q}x_{i}$ and
$\int\operatorname{d}y_{a}\mapsto\int\operatorname{d}_{t^{-1}}y_{a}$, which
would satisfy the same set of $q$-Virasoro constraints.

###  3.1 q-Virasoro constraints

In order to derive the appropriate $q$-Virasoro constraints for the rABJ model
we need to find a generalization of the difference operator $\mathcal{D}_{m}$
in (2.21) which acts both on the $x_{i}$ variables and on the $y_{a}$
variables. The fact that these two sets of variables have different couplings
$q$ and $t^{-1}$, respectively, suggests that their shifts should also be
different. By using the symmetry under the exchange of $x_{i}\leftrightarrow
y_{a}$ and $q\leftrightarrow t^{-1}$ we are lead to define the operators

$\displaystyle\mathcal{D}_{m}f(x,y)=$
$\displaystyle\sum_{i=1}^{N}D_{q,x_{i}}\left[q^{-\frac{m}{2}}x_{i}^{m+1}A_{i}(x,y;q^{-1},t^{-1})f(x,y)\right]$
(3.8)
$\displaystyle+\sum_{a=1}^{M}D_{t^{-1},y_{a}}\left[t^{\frac{m}{2}}y_{a}^{m+1}B_{a}(x,y;q^{-1},t^{-1})f(x,y)\right]$

where

$A_{i}(x,y;q,t)=\prod_{j\neq
i}\frac{tx_{i}-x_{j}}{x_{i}-x_{j}}\prod_{a=1}^{M}\frac{y_{a}-\sqrt{t/q}x_{i}}{y_{a}-q\sqrt{t/q}x_{i}}$
(3.9) $B_{a}(x,y;q,t)=\prod_{b\neq
a}\frac{q^{-1}y_{a}-y_{b}}{y_{a}-y_{b}}\prod_{i=1}^{N}\frac{x_{i}-\sqrt{t/q}y_{a}}{x_{i}-t^{-1}\sqrt{t/q}y_{a}}$
(3.10)

We now show that this are the correct difference operators to derive the rABJ
generalization of the $q$-Virasoro constraints. First we observe that the
functions $A_{i}(x,y;q,t)$ and $B_{a}(x,y;q,t)$ satisfy

$\Delta_{q,t}(x,y)^{-1}T_{q,x_{i}}\Delta_{q,t}(x,y)=q^{M}t^{-(N-1)}A_{i}(x,y;q,t)\left[T_{q,x_{i}}A_{i}(x,y;q^{-1},t^{-1})\right]^{-1}$
(3.11)

and

$\Delta_{q,t}(x,y)^{-1}T_{t^{-1},y_{a}}\Delta_{q,t}(x,y)=t^{-N}q^{M-1}B_{a}(x,y;q,t)\left[T_{t^{-1},y_{a}}B_{a}(x,y;q^{-1},t^{-1})\right]^{-1}$
(3.12)

With some straightforward algebraic manipulations we can rewrite the
constraints as

$\left\langle\mathbb{D}_{m}^{V}\Pi(x,y,{\bf t};q,t)\right\rangle=0$ (3.13)

where

$\displaystyle\mathbb{D}_{m}^{V}=$
$\displaystyle\frac{1}{1-q}\Bigg{[}\sum_{i=1}^{N}(q^{-\frac{1}{2}}x_{i})^{m}A_{i}(x,y;q^{-1},t^{-1})$
(3.14)
$\displaystyle\quad\quad\quad-q^{M+1}t^{-(N-1)}\sum_{i=1}^{N}(q^{\frac{1}{2}}x_{i})^{m}\operatorname{e}^{V(x,y)-T_{q,x_{i}}V(x,y)}A_{i}(x,y;q,t)T_{q,x_{i}}\Bigg{]}$
$\displaystyle+\frac{1}{1-t^{-1}}\Bigg{[}\sum_{a=1}^{M}(t^{\frac{1}{2}}y_{a})^{m}B_{a}(x,y;q^{-1},t^{-1})$
$\displaystyle\quad\quad\quad\quad\quad-t^{-(N+1)}q^{M-1}\sum_{a=1}^{M}(t^{-\frac{1}{2}}y_{a})^{m}\operatorname{e}^{V(x,y)-T_{t^{-1},y_{a}}V(x,y)}B_{a}(x,y;q,t)T_{t^{-1},y_{a}}\Bigg{]}$

which reduces to (2.25) when $M=0$. The last step of the derivation consists
in computing the formal adjoint of the operators $\mathbb{D}_{m}^{V}$ using
(3.5), and we find

$\displaystyle(\mathbb{D}^{V}_{m})^{\perp}=\oint_{|z|=1}\frac{\operatorname{d}z}{2\pi\operatorname{i}z}z^{-m}\frac{1}{(1-q)(1-t^{-1})}\Bigg{\\{}1-q^{M}t^{-N}\exp\left(\sum_{n=1}^{\infty}\frac{(zq^{-\frac{1}{2}})^{n}}{n}(1-t^{n})t_{n}^{\perp}\right)$
(3.15)
$\displaystyle+prz\left[1-q^{-M}t^{N}\exp\left(-\sum_{n=1}^{\infty}\frac{(q^{\frac{1}{2}}z)^{-n}}{n}(1-t^{n})t_{n}\right)\exp\left(\sum_{n=1}^{\infty}\frac{(q^{\frac{1}{2}}z)^{n}}{n}(1-t^{-n})t_{n}^{\perp}\right)\right]\Bigg{\\}}$

so that the constraints become

$(\mathbb{D}_{m}^{V})^{\perp}Z^{\mathrm{rABJ}}_{N,M}({\bf
t};q,t)=0\leavevmode\nobreak\ ,\quad\quad\quad m\geq 1$ (3.16)

A couple of observations are in order. First we remark that the operator in
(3.15) can formally be obtained from (2.28) via the substitution

$N\mapsto N_{\mathrm{eff}}:=N-\tfrac{1}{\beta}M$ (3.17)

so that we can use the $q$-Virasoro constraints to argue that there exists an
equivalence between rABJ theory for the super-group $U(N,M)$ and rCS for
$U(N-\frac{1}{\beta}M)$. The precise statement of this equivalence is given by
the identity

$\frac{Z_{N,M}^{\mathrm{rABJ}}({\bf
t};q,t)}{Z_{N,M}^{\mathrm{rABJ}}(0;q,t)}=\frac{Z_{N-\frac{1}{\beta}M}^{\mathrm{rCS}}({\bf
t};q,t)}{Z_{N-\frac{1}{\beta}M}^{\mathrm{rCS}}(0;q,t)}$ (3.18)

between normalized rCS and rABJ generating functions. Because the $q$-Virasoro
constraints in rCS are completely solvable (up to overall normalization), it
follows that the rABJ generating function is also uniquely defined by
$q$-Virasoro.

As a second observation we can ask the following question: given some specific
values of $\beta$ and $N_{\mathrm{eff}}$, are there two integers $N,M$ such
that $N_{\mathrm{eff}}=N-\frac{1}{\beta}M$? If this is the case, then we can
use the correspondence in (3.18) to define rCS theory for non-integer rank
$N_{\mathrm{eff}}$ (which does not make sense as an honest matrix model) by
instead computing the generating function of rABJ theory with integer ranks
$N,M$ as in the r.h.s. of (3.17).

We will now provide a full answer to this question in the case that
$\beta\in\mathbb{Q}$. Let $\beta=-\frac{\epsilon_{2}}{\epsilon_{1}}$ where
$\epsilon_{1},\epsilon_{2}\in\mathbb{Z}$ and we can assume without loss of
generality that $\gcd(\epsilon_{1},\epsilon_{2})=1$. We want to find
$N,M\in\mathbb{Z}$ such that

$N_{\mathrm{eff}}=N+\frac{\epsilon_{1}}{\epsilon_{2}}M$ (3.19)

If $N_{\mathrm{eff}}\in\mathbb{Z}$ we can just take $N=N_{\mathrm{eff}}$ and
$M=0$. If $N_{\mathrm{eff}}$ is not integer then the previous equation implies
that it must be a rational number. We can then assume that
$N_{\mathrm{eff}}=\frac{a}{b}$ for some coprime integers $a,b$. With these
assumptions the equation we need to solve becomes

$b\epsilon_{2}N+b\epsilon_{1}M=a\epsilon_{2}$ (3.20)

This linear Diophantine equation admits an integer solution iff there exists
an integer $e$ such that

$eb=a\epsilon_{2}\quad\quad\Leftrightarrow\quad\quad
N_{\mathrm{eff}}=\frac{a}{b}=\frac{e}{\epsilon_{2}}$ (3.21)

If $N_{\mathrm{eff}}$ has this specific form, then we can use Bézout’s lemma
to find two integers $n,m$ that solve the equation

$\epsilon_{2}n+\epsilon_{1}m=1$ (3.22)

It then follows that

$N_{\mathrm{eff}}=\frac{e}{\epsilon_{2}}=e\frac{\epsilon_{2}n+\epsilon_{1}m}{\epsilon_{2}}=(en)-\frac{1}{\beta}(em)$
(3.23)

therefore we can take $N=en$ and $M=em$ as integer solutions to the problem.

As a final remark, we observe that the choice of rational $\beta$ has a
physical interpretation in terms of Virasoro minimal models. More precisely,
the Virasoro central charge of a minimal model is parametrized by two integers
$\epsilon_{1},\epsilon_{2}$ as

$c_{\epsilon_{1},\epsilon_{2}}=1+6\frac{(\epsilon_{1}+\epsilon_{2})^{2}}{\epsilon_{1}\epsilon_{2}}=1-6\left(\sqrt{\beta}-\frac{1}{\sqrt{\beta}}\right)^{2}$
(3.24)

This suggests that there might be a connection between rABJ theories at
rational $\beta$ and Virasoro minimal models, or rather their
$q$-deformations, which could be explained by the BPS/CFT correspondence.

###  3.2 Langlands duality

The “quantum $q$-geometric” Langlands duality observed for the rCS matrix
model generalizes to the case of refined ABJ, where the action of the
involution ${{}^{\mathrm{L}}\omega}$ is much more obvious in that it now
exchanges the degrees of freedom $\\{x_{i}\\}$ with the $\\{y_{a}\\}$ while
also exchanging $q\leftrightarrow t^{-1}$ as well as $N\leftrightarrow M$.
More specifically, by using

$\displaystyle{{}^{\mathrm{L}}\omega}(t_{n})$
$\displaystyle=-\frac{1-t^{n}}{1-q^{n}}p^{\frac{n}{2}}t_{n}$ (3.25)
$\displaystyle{{}^{\mathrm{L}}\omega}(p_{n}(x,y))$
$\displaystyle=-\frac{1-t^{n}}{1-q^{n}}p^{\frac{n}{2}}p_{n}(x,y)$

we observe that not just the operators $\mathbb{D}^{V}_{m}$ and its adjoint
are manifestly invariant under ${{}^{\mathrm{L}}\omega}$ but also the actual
generating function in (3.6),

${{}^{\mathrm{L}}\omega}\left({Z^{\mathrm{rABJ}}_{N,M}({\bf
t};q,t)}\right)={Z^{\mathrm{rABJ}}_{M,N}({{}^{\mathrm{L}}\omega}({\bf
t});t^{-1},q^{-1})}={Z^{\mathrm{rABJ}}_{N,M}({\bf t};q,t)}$ (3.26)

where the first equality follows from the definition of the involution
${{}^{\mathrm{L}}\omega}$ while the second equality follows from the
definition of the generating function ${Z^{\mathrm{rABJ}}_{N,M}({\bf
t};q,t)}$. Notice that this is not the case for the generating function of
rCS, where the symmetry is not manifest at the level of the integral but it
becomes so only when the $q$-Virasoro solution is normalized by the partition
function.

Moreover, we can use this improved understanding of Langlands duality in rABJ
theory to explain the action of ${{}^{\mathrm{L}}\omega}$ in rCS. To do so, we
use the identification (3.7) to which we apply ${{}^{\mathrm{L}}\omega}$ to
obtain

${{}^{\mathrm{L}}\omega}\left(\frac{Z^{\mathrm{rCS}}_{N}({\bf
t};q,t)}{Z^{\mathrm{rCS}}_{N}(0;q,t)}\right)={{}^{\mathrm{L}}\omega}\left(\frac{Z^{\mathrm{rABJ}}_{N,0}({\bf
t};q,t)}{Z^{\mathrm{rABJ}}_{N,0}(0;q,t)}\right)=\frac{Z^{\mathrm{rABJ}}_{0,N}({{}^{\mathrm{L}}\omega}({\bf
t});t^{-1},q^{-1})}{Z^{\mathrm{rABJ}}_{0,N}(0;t^{-1},q^{-1})}=\frac{Z^{\mathrm{rCS}}_{-\beta
N}({{}^{\mathrm{L}}\omega}({\bf t});t^{-1},q^{-1})}{Z^{\mathrm{rCS}}_{-\beta
N}(0;t^{-1},q^{-1})}$ (3.27)

where in the last equality we used (3.18). This provides a better explanation
for the transformation rule of the rCS matrix model.

###  3.3 Averages of Super-Macdonald polynomials

Similarly to the case of rCS, we can study superintegrability of characters
for rABJ theory. If we regard rABJ as the super-group version of rCS we can in
principle consider expectation values of any irreducible character of
$U(N|M)$. These correspond to insertions of unknot Wilson loop in irreducible
representations of the super-group and these were already considered in [40,
42] for the unrefined case. There it was argued that on $S^{3}$ the $U(N|M)$
theory can be Higgsed down to $U(N-M)$ so that the expectation value of a
Wilson loop labeled by a maximally atypical representation of the super-group
is equal to the expectation value of the corresponding Wilson loop for the
bosonic theory at the effective rank $N-M$. We argue that this symmetry
breaking phenomenon should admit a refinement which is the physical
manifestation of the correspondence in (3.18).

In this paper we do not consider arbitrary irreducible representations of
$U(N|M)$, but we restrict to those whose characters can be written as
polynomial combinations of the deformed power-sums in (3.3). These characters
belong to the ring of symmetric function in the variables $\\{x_{i}\\}$ and
$\\{y_{a}\\}$, which is known to be generated by a basis of super-Macdonald
polynomials $\operatorname{SMcdP}_{\lambda}(x,y)$ defined in [41, 38]. These
polynomials are obtained from the usual Macdonald functions by specializing
the power-sums $p_{n}$ to the deformed power-sums $p_{n}(x,y)$, i.e.

$\operatorname{SMcdP}_{\lambda}(x,y)=\operatorname{McdP}_{\lambda}\left(\\{p_{n}=p_{n}(x,y)\\}\right)$
(3.28)

Now, using that

$\langle p_{n}(x,y)\rangle=t^{\perp}_{n}Z^{\mathrm{rABJ}}_{N,M}({\bf
t},q,t)|_{{\bf t}=0}$ (3.29)

we can immediately generalize the superintegrability formula (2.32) by
substituting the rank $N$ with the effective rank $N-\tfrac{1}{\beta}M$ to
give

$\frac{\langle\operatorname{SMcdP}_{\lambda}(x,y)\rangle}{\langle
1\rangle}=\frac{\operatorname{McdP}_{\lambda}\left(\left\\{p_{n}=\frac{-(-pq^{\frac{1}{2}}rt^{N}q^{-M})^{n}}{1-t^{n}}\right\\}\right)}{\operatorname{McdP}_{\lambda}\left(\left\\{p_{n}=\frac{1}{1-t^{n}}\right\\}\right)}{\operatorname{McdP}_{\lambda}\left(\left\\{p_{n}=\frac{1-t^{nN}q^{-nM}}{1-t^{n}}\right\\}\right)}$
(3.30)

The rABJ generating function admits a similar character expansion to that in
(2.39), therefore knowledge of the expectation value of Super-Macdonald
polynomials implies full knowledge of the function
$Z^{\mathrm{rABJ}}_{N,M}({\bf t};q,t)$.

##  4 Conclusions and discussions

In this paper we considered matrix models for a refinement of CS theory and
ABJ theory and we showed that there is an action of the $q$-Virasoro algebra
on their generating functions of observables. Moreover, we argued that the
condition that such generating functions are annihilated by all positive
$q$-Virasoro generators, provides an infinite set of homogeneous constraints
which admit a unique solution in the space of formal power series in higher
times. We derived such solutions both by recursion on the correlation
functions and also as character expansion formulas in terms of Macdonald
averages. Our methods provide an alternative point of view on the results of
[2, 11] for rCS while they give new interesting predictions for the less-
studied refinement of ABJ theory (and ABJM when $N=M$).

Here we list some related questions and future research directions.

  * •

The first and perhaps most obvious question one can ask is whether the results
of this paper can be extended to the case of CS theory in the presence of
matter fields. In the series of papers [14, 43, 8] the case of 3d
$\mathcal{N}=2$ Yang-Mills with an adjoint matter field and $N_{f}$
fundamental (anti)-chirals was considered and $q$-Virasoro constraints were
derived and solved. However, a strong technical requirement for the complete
solvability of the equations was that the effective CS level be set to zero.
In the present paper the situation is somewhat reversed in the sense that the
level can be chosen arbitrarily but there are no matter fields. It is not
clear to us if these are just some technical difficulties or whether there is
a deeper relation between uniqueness of the solution of the constraints and
the specific choice of parameters and field content of the 3d gauge theory.

  * •

Another natural extension of our results would be that of considering
$q$-Virasoro constraints in matrix models for more general knots, such as the
torus knot matrix models. Even though the approach of [2, 11] using the
explicit form of the Kac-Peterson modular matrices has been successful in
computing refined knot invariants for torus knots, it is still not clear how
to write down a refined matrix model in that case. We leave the investigation
of $q$-Virasoro for torus knots matrix models for future projects, however we
do make some remarks of the relation between torus knot invariants and unknot
invariants in Appendix B, alas, for the unrefined case.

  * •

In [44] an elliptic deformation of the fusion ring of $\widehat{su}(N)_{k}$
was introduced and the corresponding elliptic deformations of the modular $S$
and $T$ matrices were also computed. This deformation is known to reproduce
trigonometric/Macdonald deformation when the elliptic parameter is set to
zero. It would be interesting to investigate whether there exists an elliptic
Virasoro [45] version of the constraints for the corresponding elliptic
deformation of CS matrix model and whether those constraints are uniquely
solvable as in the case of the present paper.

  * •

Finally, a central question raised by our results is whether there exists a
wider class of potentials $V(x)$ such that the $q$-Virasoro constraints can be
solved fully and uniquely in terms of correlators or averages of characters.
Namely, in [14, 8] it was shown that potentials of the form
$V(x)=\sum_{k=1}^{N_{f}}\operatorname{Li}_{2}(u_{k}x;q)$ satisfying

$\operatorname{e}^{V(x)-V(qx)}=\prod_{k=1}^{N_{f}}(1-u_{k}x)^{-1}$

lead to a full solution of the constraints for $N_{f}=1,2$, where $N_{f}$ is
the step of the recursion. For higher $N_{f}$ it was shown that the
constraints cannot fix all correlators but only a subset and that the
remaining correlators do depend on additional non-canonical choices of initial
data for the recursion. In this paper we instead argue that there is another
potential which also leads to constraints that have a unique solution only
depending on a choice of normalization. The $q$-shift transformation of the
rCS potential $V(x)=\tfrac{\log^{2}x}{2\log q}$ is

$\operatorname{e}^{V(x)-V(qx)}=q^{-\frac{1}{2}}x^{-1}$

and it induces a recursion of step 1. It appears from this discussion that
whenever

$\operatorname{e}^{V(x)-V(qx)}=(\text{polynomial in $x$ of degree $d$})^{-1}$

then the $q$-Virasoro constraints can be written as recursion relations
between correlators of order $m$ and correlators of order $m-d$, which appear
to be uniquely solvable only for $d=1$ and $d=2$. It is not clear to us
whether these are the only cases for which the solution of the constraints is
unique.

## Acknowledgments

All authors are supported in part by the grant “Geometry and Physics” from the
Knut and Alice Wallenberg foundation.

##  Appendix A Symmetric functions

For completeness we collect some basic notions on the theory of symmetric
functions. See [17] for more details.

We start by recalling some definitions about partitions. An integer partition
$\lambda$ is a sequence of non-increasing integer numbers

$\lambda=[\lambda_{1},\lambda_{2},\dots]$ (A.1)

with a finite number of non-zero terms. We denote by $\ell(\lambda)$ the
length of the partition, defined as

$\ell(\lambda)=\mathrm{Card}\\{j:\lambda_{j}\neq 0\\}$ (A.2)

and we denote $|\lambda|$ the weight or size or degree of the partition,

$|\lambda|=\sum_{i\geq 1}\lambda_{i}$ (A.3)

and we say that $\lambda\vdash n$ if $|\lambda|=n$. The conjugate partition to
$\lambda$ is denoted by $\lambda^{\prime}$. For later convenience we also
introduce the staircase partition ${\boldsymbol{\delta}}$ defined as

${\boldsymbol{\delta}}=[N-1,N-2,\dots,1,0]$ (A.4)

We define $\operatorname{Aut}(\lambda)$ as the group of automorphisms of a
given partition, and its cardinality can be computed as

$|\operatorname{Aut}(\lambda)|=\prod_{i\geq
1}\mathrm{Card}\\{j:\lambda_{j}=i\\}!$ (A.5)

It is also useful to define the function

$z_{\lambda}=|\operatorname{Aut}(\lambda)|\prod_{i\geq 1}\lambda_{i}$ (A.6)

Additionally we define a scalar product on partitions as

$\lambda\cdot\mu=\sum_{i\geq 1}\lambda_{i}\,\mu_{i}$ (A.7)

Symmetric functions are elements of the ring $\Lambda$ defined as the inverse
limit of the graded ring $\Lambda_{N}=\mathbb{Q}[x_{1},\dots,x_{N}]^{S_{N}}$
of symmetric polynomials in $N$ variables with respect to the projections
$\Lambda_{M}\to\Lambda_{N}$ for $M\geq N$. The ring $\Lambda$ has many bases,
among which the elementary symmetric functions $e_{n}$, the complete
homogeneous symmetric functions $h_{n}$, the power-sums symmetric functions
$p_{n}$ and the Schur symmetric functions $\operatorname{Schur}_{\lambda}$.
Upon defining the power-sums as

$p_{n}=\sum_{i\geq 1}x_{i}^{n}\leavevmode\nobreak\ ,\quad\quad\quad
p_{\lambda}=\prod_{i\geq 1}p_{\lambda_{i}}$ (A.8)

we can obtain all other basis via polynomial combinations, e.g.

$h_{m}=\sum_{\lambda\vdash
m}{z_{\lambda}^{-1}}{p_{\lambda}}\leavevmode\nobreak\ ,\quad\quad\quad
e_{m}=\sum_{\lambda\vdash
m}(-1)^{\ell(\lambda)}{z_{\lambda}^{-1}}{p_{\lambda}}$ (A.9)
$\operatorname{Schur}_{\lambda}=\det_{1\leq
i,j\leq\ell(\lambda)}\left[h_{\lambda_{i}-i+j}\right]$ (A.10)

Moreover, there exists an inner product on symmetric functions defined by the
relations

$(p_{\lambda},p_{\mu})=z_{\lambda}\delta_{\lambda\mu}$ (A.11)

or equivalently

$(\operatorname{Schur}_{\lambda},\operatorname{Schur}_{\mu})=\delta_{\lambda\mu}$
(A.12)

Let $\Lambda_{q,t}=\Lambda\otimes\mathbb{F}_{q,t}$ with $\mathbb{F}_{q,t}$ the
field of rational functions in $q$ and $t$. Then there exists a deformed inner
product defined by

$(p_{\lambda},p_{\mu})_{q,t}=z_{\lambda}(q,t)\delta_{\lambda\mu}\leavevmode\nobreak\
,\quad\quad\quad z_{\lambda}(q,t)=z_{\lambda}\prod_{i\geq
1}\frac{1-q^{\lambda_{i}}}{1-t^{\lambda_{i}}}$ (A.13)

Macdonald functions are defined as the $q,t$-deformation of Schur symmetric
functions which are orthogonal w.r.t. the inner product $(\cdot,\cdot)_{q,t}$.
More precisely, there are two dual basis of Macdonald functions,
$\operatorname{McdP}_{\lambda}$ and $\operatorname{McdQ}_{\mu}$ with

$(\operatorname{McdP}_{\lambda},\operatorname{McdQ}_{\mu})_{q,t}=\delta_{\lambda\mu}$
(A.14)

For $f\in\Lambda_{q,t}$ we use the plethystic notation $f(\\{p_{n}\\})$ to
indicate that it is a symmetric function written as a polynomial in the power-
sums, while we use the notation $f(x)\equiv f(x_{1},x_{2},\dots)$ to indicate
that it is an actual symmetric polynomial in the variables $\\{x_{i}\\}$.

##  Appendix B Remarks on torus knots

In the main part of the paper we considered matrix models for unknot Wilson
loops. In this section we provide some remarks for the case of torus knots.
This is the only other class of knots for which a matrix model description has
been derived [46, 47, 48, 49], unfortunately no refinement of this matrix
model has been worked out yet. The general expectation is that there should be
a matrix model computing knot polynomials for any kind of knot, and perhaps
each of such matrix model should admit a refinement.

It would be interesting to ask whether there exists a $q$-Virasoro action on
the generating function of Wilson loop expectation values in all of these
matrix models666See [50] for earlier attempts of using classical Virasoro
constraints to compute correlators in unrefined torus knot matrix models.. For
the case of torus knots matrix models the $q$-difference operator in (2.21) at
($t=q$) does not seem to give rise to recursion equations that can be solved
uniquely. A more subtle modification of that operator is needed and we leave
that for future investigations. Here we make some observations on how to
employ the solution to $q$-Virasoro constraints derived in section 2.3 to make
predictions about torus knots expectation values.

The $(P,Q)$-torus knot generating function is

$\displaystyle Z^{(P,Q)}_{N}({\bf
t};q)=\frac{(PQ)^{N}}{N!}\int_{\mathbb{R}^{N}_{+}}\Delta(x^{P})\Delta(x^{Q})\prod_{i=1}^{N}\exp\left[-PQ\frac{\log^{2}x_{i}}{2\log{q}}+\nu\log
x_{i}+\sum_{n=1}^{\infty}\frac{x_{i}^{nPQ}t_{n}}{n}\right]\mathrm{d}x_{i}$
(B.1)

where $P,Q$ are positive integers defining the knot and

$\Delta(x)=\prod_{i<j}(x_{i}-x_{j})=\det_{1\leq i,j\leq
N}[x_{i}^{{\boldsymbol{\delta}}_{j}}]$ (B.2)

is the Vandermonde determinant (with ${\boldsymbol{\delta}}$ as in (A.4)).
Using the identity

$\frac{\Delta(x^{Y})}{\Delta(x)}=\frac{\det_{1\leq i,j\leq
N}[x_{i}^{{\boldsymbol{\delta}}_{j}+(Y-1){\boldsymbol{\delta}}_{j}}]}{\det_{1\leq
i,j\leq
N}[x_{i}^{{\boldsymbol{\delta}}_{j}}]}=\operatorname{Schur}_{(Y-1){\boldsymbol{\delta}}}(x)\leavevmode\nobreak\
,\quad\quad\quad Y=P,Q$ (B.3)

which follows from Weyl’s character formula for Schur polynomials, we can
write

$\displaystyle Z^{(P,Q)}_{N}({\bf t};q)=$
$\displaystyle\frac{(PQ)^{N}}{N!}\int_{\mathbb{R}^{N}_{+}}\prod_{i=1}^{N}\mathrm{d}x_{i}\,|\Delta(x)|^{2}\operatorname{Schur}_{(P-1){\boldsymbol{\delta}}}(x)\operatorname{Schur}_{(Q-1){\boldsymbol{\delta}}}(x)\times$
(B.4)
$\displaystyle\times\prod_{i=1}^{N}\exp\left[-PQ\frac{\log^{2}x_{i}}{2\log{q}}+\nu\log
x_{i}+\sum_{n=1}^{\infty}\frac{x_{i}^{nPQ}t_{n}}{n}\right]$ $\displaystyle=$
$\displaystyle(PQ)^{N}\operatorname{Schur}_{(P-1){\boldsymbol{\delta}}}(\\{u_{n}^{\perp}\\})\operatorname{Schur}_{(Q-1){\boldsymbol{\delta}}}(\\{u_{n}^{\perp}\\})\left.Z_{N}^{\mathrm{CS}}({\bf
u};q^{\frac{1}{PQ}})\right|_{u_{n}=PQ\delta_{n|PQ}t_{n/PQ}}$

where in the last line we rewrote the torus knot generating function as the
action of two adjoint Schur polynomials on the generating function of
(unrefined) CS matrix model at the coupling $q^{\frac{1}{PQ}}$. The times of
the two generating functions are identified using the fact that

$\sum_{n=1}^{\infty}\frac{x_{i}^{nPQ}t_{n}}{n}=\sum_{n=1}^{\infty}\frac{x_{i}^{n}}{n}(PQ\delta_{n|PQ}t_{{n}/{PQ}})$
(B.5)

where $\delta_{n|PQ}$ is 1 if $n$ is a multiple of $PQ$ and 0 otherwise.

Since $q$-Virasoro constraints give a full solution for
$Z_{N}^{\mathrm{CS}}({\bf u};q^{\frac{1}{PQ}})$, then we also have a
combinatorially explicit way to compute $Z^{(P,Q)}_{N}({\bf t};q)$ just by
acting with an appropriate differential operator of finite degree in times.

While it is tempting to conjecture the form of refined CS matrix model for
torus knots by deforming Schur polynomials in (B.4) to Macdonald polynomials,
we must observe that the identity (B.3) does not generalize to the Macdonald
level and therefore some other techniques are required in order to treat the
refined model.

## References

  * [1] E. Witten, “Quantum Field Theory and the Jones Polynomial,” Commun. Math. Phys. 121 (1989) 351–399.
  * [2] M. Aganagic and S. Shakirov, “Knot Homology and Refined Chern-Simons Index,” Commun. Math. Phys. 333 no. 1, (2015) 187–228, arXiv:1105.5117 [hep-th].
  * [3] N. M. Dunfield, S. Gukov, and J. Rasmussen, “The Superpolynomial for knot homologies,” arXiv:math/0505662.
  * [4] M. Khovanov, “A categorification of the Jones polynomial,” Duke Mathematical Journal 101 no. 3, (2000) 359–426, arXiv:math/9908171 [math.QA].
  * [5] M. Khovanov and L. Rozansky, “Matrix factorizations and link homology II,” Geometry & Topology 12 no. 3, (2008) 1387–1425, arXiv:math/0505056 [math.QA].
  * [6] J. Shiraishi, H. Kubo, H. Awata, and S. Odake, “A Quantum deformation of the Virasoro algebra and the Macdonald symmetric functions,” Lett. Math. Phys. 38 (1996) 33–51, arXiv:q-alg/9507034.
  * [7] A. Nedelin and M. Zabzine, “q-Virasoro constraints in matrix models,” JHEP 03 (2017) 098, arXiv:1511.03471 [hep-th].
  * [8] L. Cassia, R. Lodin, and M. Zabzine, “On matrix models and their $q$-deformations,” JHEP 10 (2020) 126, arXiv:2007.10354 [hep-th].
  * [9] O. Aharony, O. Bergman, D. L. Jafferis, and J. Maldacena, “N=6 superconformal Chern-Simons-matter theories, M2-branes and their gravity duals,” JHEP 10 (2008) 091, arXiv:0806.1218 [hep-th].
  * [10] O. Aharony, O. Bergman, and D. L. Jafferis, “Fractional M2-branes,” JHEP 11 (2008) 043, arXiv:0807.4924 [hep-th].
  * [11] M. Aganagic and S. Shakirov, “Refined Chern-Simons Theory and Knot Homology,” Proc. Symp. Pure Math. 85 (2012) 3–32, arXiv:1202.2489 [hep-th].
  * [12] M. Marino, “Chern-Simons theory, matrix integrals, and perturbative three manifold invariants,” Commun. Math. Phys. 253 (2004) 25–49, arXiv:hep-th/0207096.
  * [13] M. Aganagic, A. Klemm, M. Marino, and C. Vafa, “Matrix model as a mirror of Chern-Simons theory,” JHEP 02 (2004) 010, arXiv:hep-th/0211098.
  * [14] R. Lodin, A. Popolitov, S. Shakirov, and M. Zabzine, “Solving q-Virasoro constraints,” Lett. Math. Phys. 110 no. 1, (2020) 179–210, arXiv:1810.00761 [hep-th].
  * [15] A. J. Kirillov, “On inner product in modular tensor categories. I,” J. Amer. Math. Soc. 9 (1996) 1135–1169, arXiv:q-alg/9508017 [math.QA].
  * [16] I. Cherednik, “Macdonald’s evaluation conjectures and difference fourier transform,” Inventiones Mathematicae 122 no. 1, (Dec, 1995) 119–145, arXiv:q-alg/9412016 [math.QA].
  * [17] I. G. Macdonald, Symmetric Functions and Hall Polynomials.  Oxford Mathematical Monographs. Oxford University Press, USA, 2 ed., 1995\.
  * [18] A. Mironov and A. Morozov, “On the complete perturbative solution of one-matrix models,” Phys. Lett. B 771 (2017) 503–507, arXiv:1705.00976 [hep-th].
  * [19] A. Mironov and A. Morozov, “Correlators in tensor models from character calculus,” Phys. Lett. B 774 (2017) 210–216, arXiv:1706.03667 [hep-th].
  * [20] A. Mironov and A. Morozov, “Sum rules for characters from character-preservation property of matrix models,” JHEP 08 (2018) 163, arXiv:1807.02409 [hep-th].
  * [21] S. M. Natanzon and A. Y. Orlov, “Hurwitz numbers and BKP hierarchy,” arXiv:1407.8323 [nlin.SI].
  * [22] A. Morozov, A. Popolitov, and S. Shakirov, “On (q,t)-deformation of Gaussian matrix model,” Phys. Lett. B 784 (2018) 342–344, arXiv:1803.11401 [hep-th].
  * [23] E. Frenkel and N. Reshetikhin, “Deformations of $\mathcal{W}$-algebras associated to simple Lie algebras,” Comm. Math. Phys. 197 (1998) 1–32, arXiv:q-alg/9708006 [math.QA].
  * [24] E. Frenkel and D. Hernandez, “Langlands duality for finite-dimensional representations of quantum affine algebras,” Lett. Math. Phys. 96 (2011) 217–261, arXiv:0902.0447 [math.QA].
  * [25] M. Aganagic, E. Frenkel, and A. Okounkov, “Quantum $q$-Langlands Correspondence,” Trans. Moscow Math. Soc. 79 (2018) 1–83, arXiv:1701.03146 [hep-th].
  * [26] P. Bouwknegt and K. Pilch, “The Deformed Virasoro Algebra at Roots of Unity,” Comm. Math. Phys. 196 (1998) 249–288, arXiv:q-alg/9710026 [math.QA].
  * [27] A. Nigro, “A q-Virasoro Algebra at roots of unity, Free Fermions and Temperley Lieb Hamiltonians,” J. Math. Phys. 57 no. 4, (2016) 041702, arXiv:1211.1067 [math-ph].
  * [28] H. Itoyama, T. Oota, and R. Yoshioka, “2d-4d Connection between $q$-Virasoro/W Block at Root of Unity Limit and Instanton Partition Function on ALE Space,” Nucl. Phys. B 877 (2013) 506–537, arXiv:1308.2068 [hep-th].
  * [29] H. Itoyama, T. Oota, and R. Yoshioka, “$q$-Virasoro/W Algebra at Root of Unity and Parafermions,” Nucl. Phys. B 889 (2014) 25–35, arXiv:1408.4216 [hep-th].
  * [30] A. Kapustin and B. Willett, “Wilson loops in supersymmetric Chern-Simons-matter theories and duality,” arXiv:1302.2164 [hep-th].
  * [31] A. Morozov, A. Popolitov, and S. Shakirov, “Quantization of Harer-Zagier formulas,” Phys. Lett. B 811 (2020) 135932, arXiv:2008.09577 [hep-th].
  * [32] A. Morozov, A. Popolitov, and S. Shakirov, “Harer-Zagier formulas for knot matrix models,” Phys. Lett. B 818 (2021) 136370, arXiv:2102.11187 [hep-th].
  * [33] J. Harer and D. Zagier, “The Euler characteristic of the moduli space of curves,” Inventiones mathematicae 85 (1986) 457–485.
  * [34] P. J. Forrester, “Global and local scaling limits for the $\beta=2$ Stieltjes–Wigert random matrix ensemble,” arXiv:2011.11783 [math-ph].
  * [35] N. Drukker and D. Trancanelli, “A Supermatrix model for N=6 super Chern-Simons-matter theory,” JHEP 02 (2010) 058, arXiv:0912.3006 [hep-th].
  * [36] B. Eynard and T. Kimura, “Toward $\mathrm{U}(N|M)$ knot invariant from ABJM theory,” Lett. Math. Phys. 107 no. 6, (2017) 1027–1063, arXiv:1408.0010 [hep-th].
  * [37] F. Nieri, Y. Pan, and M. Zabzine, “3d Expansions of 5d Instanton Partition Functions,” JHEP 04 (2018) 092, arXiv:1711.06150 [hep-th].
  * [38] F. Atai, M. Hallnäs, and E. Langmann, “Super-Macdonald polynomials: Orthogonality and Hilbert space interpretation,” arXiv:2103.07400 [math.QA].
  * [39] T. Kimura and F. Nieri, “Intersecting Defects and Supergroup Gauge Theory,” arXiv:2105.02776 [hep-th].
  * [40] V. Mikhaylov and E. Witten, “Branes And Supergroups,” Commun. Math. Phys. 340 no. 2, (2015) 699–832, arXiv:1410.1175 [hep-th].
  * [41] A. N. Sergeev and A. P. Veselov, “Deformed macdonald-ruijsenaars operators and super macdonald polynomials,” Commun. Math. Phys. 288 (2009) 653–675, arXiv:0707.3129 [math.QA].
  * [42] V. Mikhaylov, “Analytic Torsion, 3d Mirror Symmetry And Supergroup Chern-Simons Theories,” arXiv:1505.03130 [hep-th].
  * [43] L. Cassia, R. Lodin, A. Popolitov, and M. Zabzine, “Exact SUSY Wilson loops on S3 from $q$-Virasoro constraints,” JHEP 12 (2019) 121, arXiv:1909.10352 [hep-th].
  * [44] J. F. van Diejen and T. Görbe, “Elliptic Ruijsenaars difference operators, symmetric polynomials, and Wess-Zumino-Witten fusion rings,” arXiv:2106.14919 [math.QA].
  * [45] F. Nieri, “An elliptic Virasoro symmetry in 6d,” Lett. Math. Phys. 107 no. 11, (2017) 2147–2187, arXiv:1511.00574 [hep-th].
  * [46] R. Lawrence and L. Rozansky, “Witten–Reshetikhin–Turaev Invariants of Seifert Manifolds,” Communications in Mathematical Physics 205 (1999) 1432–0916.
  * [47] C. Beasley, “Localization for Wilson Loops in Chern-Simons Theory,” Adv. Theor. Math. Phys. 17 no. 1, (2013) 1–240, arXiv:0911.2687 [hep-th].
  * [48] J. Kallen, “Cohomological localization of Chern-Simons theory,” JHEP 08 (2011) 008, arXiv:1104.5353 [hep-th].
  * [49] A. Brini, B. Eynard, and M. Marino, “Torus knots and mirror symmetry,” Annales Henri Poincare 13 (2012) 1873–1910, arXiv:1105.2012 [hep-th].
  * [50] O. Dubinkin, “On the Virasoro constraints for torus knots,” J. Phys. A 47 no. 48, (2014) 485203, arXiv:1307.7909 [hep-th].