On matrix models and their $q$-deformations

Motivated by the BPS/CFT correspondence, we explore the similarities between the classical $\beta$-deformed Hermitean matrix model and the $q$-deformed matrix models associated to 3d $\mathcal{N}=2$ supersymmetric gauge theories on $D^2\times_{q}S^1$ and $S_b^3$ by matching parameters of the theories. The novel results that we obtain are the correlators for the models, together with an additional result in the classical case consisting of the $W$-algebra representation of the generating function. Furthermore, we also obtain surprisingly simple expressions for the expectation values of characters which generalize previously known results.


Introduction
The BPS/CFT correspondence [1,2] has provided a mapping between the exact values of partition functions and certain BPS observables for supersymmetric theories on the one hand, to conformal field theories in two dimensions on the other hand. The computation of partition functions and BPS observables has been aided through the programme called localization (as reviewed in [3]), where the evaluation of infinite dimensional integrals reduces to evaluation only at specific points. Such localization computations of partition functions sometimes result in the form of a finite dimensional matrix model, and it is such partition functions that will be of interest here.
In the light of this BPS/CFT correspondence, we derive and solve the Virasoro constraints that generating functions (τ -function) for certain classical and "quantum" models satisfy. These constraints can be considered a type of Ward identities and were first studied in [4,5]. Here, the notion of "quantum" will be related to the introduction of a special class of q-deformations with respect to an un-deformed model. In the classical case, we will be considering the β-deformed Hermitean matrix model, in other words a one parameter deformation of the standard Hermitean matrix model [6,7]. On the quantum side, we will explore the 3d N = 2 supersymmetric theory with U (N ) gauge group on both D 2 × q S 1 [8,9] and also S 3 b [10,11] with one adjoint chiral and an arbitrary number of anti-chiral fundamental multiplets. The model on D 2 × q S 1 has also been referred to as the (q, t)-model [12] where the parameters q and t are the two deformation parameters. In order to study the gauge theory on S 3 b we will use a construction which has been called the modular double [13]. The name refers to the picture of gluing two instances of D 2 × q S 1 to obtain S 3 b , something which is also mirrored in the algebraic structure at the level of partition functions. This modular double property of the S 3 b partition function will be alluded to in Section 3. In the gauge theory examples these partition functions are expressed as matrix models, originating from a localization computation which we here simply assume the result of. Then, we both derive the constraint equations that these models satisfy explicitly and we also show how the resulting constraint can be solved in a recursive fashion. In other words we illustrate how any correlator of the theory can be determined using a finite number of steps of the recursion relation. In the case of the β-deformed Hermitean matrix model, we could in addition to the correlators also find the W -algebra representation of the generating functions. This is a representation in which the generating function is expressed through the action of a single operator acting on a simple function. Thus, the results which are novel here for the classical case are the correlators (presented in (2.36) and (2.46)) and the W -representations of generating functions (in (2.35) and (2.45)). This generalize the result of [14][15][16] by introducing additional parameter dependence. In the case of the gauge theories on D 2 × q S 1 and S 3 b , the results are in terms of correlators (given in (3.54) and (3.74)), and they are extending the results of [17] by introducing another deformation parameter.
We also comment on the fact that averages of certain functions, when computed with respect to the measure of the partition function in question, take a particularly simple form. It is worth noting that this simplification is not expected a priori. These special functions are the Schur polynomials in the case of the standard Hermitean matrix model, Jack polynomials in the case of the β-deformed Hermitean matrix model and finally Macdonald polynomials in the case of the gauge theories on D 2 × q S 1 and S 3 b . The existence of such formulas has been referred to as the property of super-integrability of the model [16,18]. In the classical case we present formulas for the averages of Jack polynomi-als (in (2.37) and (2.47)), and in the quantum case we give the formulas for averages of Macdonald polynomials (in (3.69) and (3.82)) which improve and extend the results of [19].
To further clarify the relation between the classical and the two quantum models we have in mind, we can illustrate the relations between the models as shown in Figure 1. Here we show the various deformations and limits to obtain one model from the other, together with the corresponding polynomial (whose average has a simple formula) for each model. Furthermore, we can also perform a matching between the parameters of the models as follows. The parameter β of the classical Hermitean matrix model can be related to the mass t of the adjoint chiral in the quantum model. The polynomial degree p of the potential V in the classical model can be related to the number of fundamental anti-chiral fields N f in the quantum model. Then we can also match the coupling constants a k appearing in the classical potential V with the masses of the fundamental anti-chiral fields u k . Macdonald Jack Schur ×2 q-deformation Classical limit t = q β , q → 1 β-deformation Schur limit β → 1 The outline of the paper is as follows. In Section 2 we begin with reviewing the basics of the simplest matrix model, the Hermitean 1-matrix model, and then show how the Virasoro constraints (or Ward identities) for the model are derived. We then show how to solve these constraints using a recursive procedure, where some of the results generalize previously known results. In Section 3 we then move on to describe what could be considered quantum versions of the models outlined in Section 2. Moreover, these models correspond to certain supersymmetric gauge theories and in particular 3d N = 2 theories with U (N ) gauge group on D 2 × q S 1 or S 3 b . Similarly to the previous section, we derive the q-Virasoro constraints which these models satisfy, and also show how to solve the constraints recursively to obtain novel results for the correlators of the models. A semi-classical expansion is also presented in order to match with the corresponding classical matrix model. Then in Section 4, we conclude and suggest directions for further study. The details of special functions and of symmetric functions are left to Appendices A and B respectively. In Appendix C, we discuss the relation between the constraint operators and the generators of the q-Virasoro algebra and in Appendix D we perform the analysis of the asymptotic behaviour and convergence of the S 3 b partition function.

Review of the classical models
In this section we set the stage for the definition of the "quantum" matrix models by first reviewing the main features of the classical Hermitean matrix model (see [20,21] for an introduction to the topic). Here we present various well-known results by restating them in a way that makes it straightforward to match with the corresponding q-deformed case discussed in Section 3. Moreover, we extend and improve upon the previously known formulas for the W -representation of the generating function.

Definitions
Let us begin by recalling some details about classical matrix models. The simplest example of a matrix model is the Hermitean 1-matrix model (reviewed in [6]), whose degrees of freedom are represented by the Hermitean N × N matrix Φ. The observables of the theory are the traces Tr Φ s of the basic field and their expectation values can be neatly encoded into a generating function defined as Here V (Φ) is a complex function (usually a polynomial) called the potential while {τ s } are conjugate variables to the traces and are usually referred to as the time variables collectively denoted by τ . The integral is over the domain H N which is taken to be the space of all N × N Hermitean matrices while the measure d Φ is the standard Lebesgue measure on H N which is invariant under conjugation by unitary matrices. The generating function (2.1) is regarded as a formal power series in the times, however the coefficients of this expansion are integral functions and as such they must be convergent over the domain of integration and in some region of parameter space. For arbitrary potential V there can be analytical issues with defining these integrals and for any specific choice one should perform an in-depth study.
For a function O ∶ H N → C, we define its (un-normalized) expectation value or quantum average as and for convenience we also define a time-dependent expectation value by inserting O in the generating function as Derivatives of the generating function Z(τ ) w.r.t. the times then compute the timedependent expectation values of all possible single and multi-trace operators in the field Φ. Sending all the times to zero then yields the corresponding time-independent average, As a function of the times {τ s }, Z(τ ) should be regarded as formal power series which we can expand as where we rewrote the series as a sum over integer partitions ρ of arbitrary size, and we also defined the correlation functions c ρ as the expectation values of multi-trace operators whose powers are specified by the partition ρ = (ρ 1 , . . . , ρ ), namely The empty correlator c ∅ is by definition equal to the partition function Z = Z(0).
If the potential is an invariant function under conjugation of the argument by a unitary matrix, then one can use the adjoint action of U (N ) over H N to diagonalize Φ and rewrite the generating function as an integral over the eigenvalues {λ i }. The Lebesgue measure splits as the product of a Vandermonde determinant ∆(λ), the flat measure ∏ i d λ i over the space of eigenvalues and the Haar measure of U (N ), Up to a constant overall factor, we can then write the integral as It is common at this point to introduce a 1-parameter deformation of the model by generalizing the usual Vandermonde term as follows (for details we refer to [7]) where β is a positive integer number (which can be analytically continued to the complex plane). Finally, the generating function becomes (2.10) In the limit β → 1 one recovers the un-deformed model in (2.8). In the following we will always assume a β-deformation and therefore we will drop the label on the generating function.
For the purpose of our computations we are interested in the specific case of a potential function which is polynomial in the eigenvalues and takes the explicit form which depends on the integer parameter p and on the complex numbers a k which can be regarded as inverse coupling constants 1 . It is also worth mentioning that this kind of potential is special in the sense that we can obtain it via a constant shift of the time variables as τ s ↦ τ s − a s s, s = 1, . . . , p . (2.12) Assuming this potential, we will then for clarity introduce an index p on the generating function in (2.10) and the expectation value (2.3), i.e. Z p (τ ) and ⟨O⟩ p τ . Observe that while we assume that the dependence on the time variables {τ s } is only formal, after the shift (2.12), we need to carefully study the functional dependence of the generating function on the parameters a k , and in particular we need to make sure that the integral in (2.10) does indeed converge. In the eigenvalue model of (2.10) the contour of integration is taken to be the real domain R N being the range of the eigenvalues of an Hermitean matrix in H N , however when the potential V (λ i ) is introduced one must modify the contour in such a way that the integral is still convergent, possibly complexifying the variables λ i . For p ≠ 2 for instance, if Re(a p ) > 0 the integral is convergent and well-defined but only over half the real line 2 , i.e. for positive eigenvalues, while for p = 2 (and Re(a 2 ) > 0) the integral makes sense over the whole real N -dimensional space R N . In general the Ward identities do not depend on a specific choice of contour (provided there are no additional boundary terms) and one can regard different contours as different branches of the partition function, corresponding to different phases of the theory.
For the special case p = 1 we also remark that the model we described has a close relative in the complex 1-matrix model [22,23] for a complex N × N matrix M and its adjoint M † . Upon the change of variable to the Hermitean matrix Φ = M M † , we can re-write the generating function as an integral over the positive eigenvalues λ i . Taking the potential V to be a quadratic function 3 , we can write the most general form of this model as that of the (β-deformed) Wishart-Laguerre eigenvalue model (reviewed for instance in [24]) where ν is an additional parameter corresponding to the insertion of a determinant term of the form (det M M † ) ν . The integral is convergent provided the power of the determinant satisfies Re(ν) > −1 .
Formally, we can reabsorb this term inside of the potential V by writing it as a logarithmic interaction Even though such determinant insertions are degenerate in the usual Hermitian matrix model because they make the Virasoro constraints ill-defined, in the case of p = 1, as we shall show in the next section, this insertion is allowed and indeed gives an additional 1-parameter deformation which has a direct counterpart in the quantum case.
Let us pause here to make a comment on conventions and notation. Since the matrix model is built out of invariant functions w.r.t. the adjoint action of U (N ), we have that upon diagonalization and rewriting the generating function as an integral over the eigenvalues, there is a residual S N Weyl symmetry that permutes the variables {λ i }. It is a well known fact that the ring of symmetric functions has a basis given by the power-sum variables {p s } defined as which are precisely the variables that couple to the times {τ s } in the generating function. Derivatives with respect to the time τ s correspond to the insertions of p s . Another useful fact about symmetric function is that there exists an orthonormal basis provided by the Schur functions. The elements of this linear basis are symmetric polynomials labeled by integer partitions and are in 1-to-1 correspondence with the linear characters of U (N ). In the following we will often use that Schur polynomials can be expressed through powersums, for example for all partitions of degree 3 and lower. We refer to Appendix B for more details on the subject.
When one computes averages of Schur functions, these averages sometimes take a unexpectedly simple form as observed in [19]. For instance in the case of p = 2 and when a 1 = 0 we have the expression In what follows we will give examples of such averages.

Virasoro constraints
It is often very useful in QFT to consider Ward identities for the path integral of the theory at hand. In the case of matrix models the QFT is 0-dimensional and the Ward identities have a very clear differential geometric interpretation. Let the partition function be described as the integral over a domain X of the differential form Ω. If we consider an infinitesimal diffeomorphism generated by the vector field ξ over X, then we can use the vector to deform infinitesimally the form Ω and then compute the integral of the variation as the Lie derivative, namely δZ = Since Ω is a top form on X, we can write the Lie derivative as an exact form therefore its integral on X can only receive contribution by evaluating the form ι ξ Ω at the boundary of X (by Stokes theorem). Assuming that this form vanishes at the boundary, we get a non-trivial constraint equation corresponding to the fact that the variation δZ is identically zero.
In the case of the matrix models in Section 2.1 there is a natural family of vector fields given by which correspond to the generators of a Virasoro Lie algebra V ir diagonally embedded into V ir N , so that the vectors ξ n are invariant under permutations of the coordinates λ i . The differential form Ω is the integrand in (2.10) and the integration domain is X = R N . Writing explicitly the top form as Ω = f (λ) ∏ N i=1 d λ i , one can compute the total variation in (2.21) as which upon integration together with the notation in (2.3) leads to the Ward identity (2.24) These constraint equations are called the Virasoro constraints. As mentioned above, the determinant insertion depending on the parameter ν is only allowed in the case of p = 1. What one does at this point is to rewrite these expectation values as derivatives in times using the identity whenever all the powers s 1 , . . . , s k are non-negative (if some of the s l = 0 we just substitute the corresponding sum with multiplication by N ). This can be done for all n ≥ −1 if ν = 0, but for ν ≠ 0 the n = −1 constraint cannot be rewritten as a partial differential equation. The final form of the Virasoro constraints is then where U n is the differential operator that implements the n-th constraint and the operators L n are the standard generators of the Virasoro algebra 4 defined as (2.27) By construction then, (2.26) states that the generating function Z p (τ ) is in the common kernel of all such operators U n . In the following sections we study the properties of this kernel.

Solving the constraints
A legitimate question one might ask at this point is how strong are the Virasoro constraints in (2.26). Can they be used to determine the generating function Z p (τ ) and if so, what is the degeneracy of the solution? The answer to these questions was found in [14] via a W -algebra representation, which states that the solution is essentially unique if p = 1, 2 while for p ≥ 3 there is a degeneracy in the space of solutions which allows to determine Z p (τ ) only when additional information on the correlation functions is provided [25]. More recently in [26] the solution for p = 1, 2 was also found in the β-deformed model. We will now review the details of the derivation of the solution to the constraints and the issues one encounters when such a unique solution does not exists.

p = 1
The case p = 1 is the only one that admits a determinant insertion as discussed in (2.14), thus in what follows we will always assume dependence on ν in the case of p = 1. The n = −1 constraint in (2.26) is not well defined as a differential equation as it contains both negative powers of the {λ i } as well as additional boundary terms. For these reasons we restrict ourselves to consider only Virasoro constraints for n ≥ 0 where from now on we use ∂ n = ∂ ∂τn to ease notation. To obtain the solution we re-sum all constraints to construct the operator which we have rewritten as the difference of two operators: D = ∑ ∞ s=1 sτ s ∂ s is the dilatation operator and is a "shifted" W -algebra generator also called cut-and-join operator. Here the word shifted refers to the fact that W −1 is of degree 1 with respect to the grading introduced by the operator D (i.e. [D, W −1 ] = W −1 ).
Let us analyze the properties of these operators. First we remark that they are linear operators acting on the infinite dimensional vector space underlying the commutative ring This vector space has a natural basis over C given by the monomials, i.e. products of times labeled by integer partitions ρ a∈ρ τ a . (2.31) The ordering of the basis of the vector space is the one induced by the ordering on integer partitions, namely ordering by degree and lexicographic ordering between partitions of equal degree With these conventions in place one can show that U is triangular and that D is its diagonal part while W −1 is its off-diagonal part. More precisely, one finds that U is triangular also with respect to the weaker partial order induced by the monomial degree only (partition size). Moreover, we have that D acts on monomials as multiplication by the degree of the corresponding partition, D a∈ρ τ a = a∈ρ a a∈ρ τ a = ρ a∈ρ τ a (2.33) so that det U = det D = 0 because there is one zero-eigenvalue corresponding to the empty partition. However, since all other eigenvalues are non-zero, the kernel of U is exactly 1-dimensional. This means that the generating function Z p=1 (τ ), regarded as a vector in this space, is uniquely defined up to a constant multiplicative factor which corresponds to the normalization of the trivial correlator c ∅ .
The full solution can be derived by recursively solving the equation where ρ ⊢ d denotes that ρ is an integer partition of d with d being the degree in times. Then, using the fact that W −1 is of degree 1, we can write A full solution of the p = 1 model, up to degree 3 is given by the correlators (2.36) Observe that all correlators of degree higher than 1 are proportional to c {1} , which is a consequence of the fact that deg(W −1 ) = 1 and that there is only 1 partition in degree 1.
Our solution of the p = 1 model is slightly more general than the one in [15,16] as we allow for a determinant deformation of parameter ν. For the special case ν = 0 for the correlators above, we recover the formulas of [15].

Averages of characters
Another remarkable property of this model is that of super-integrability [16,18], meaning that there are some observables whose expectation values satisfy a particularly nice formula. Namely, one observes that expectation values of characters can be expressed as simple combinations of the same characters evaluated at some specific "points" (see [19] for the original observation of this fact).
In the case of the β-deformed model, the natural characters to consider are the 1-parameter family of symmetric polynomials called Jack polynomials Jack ρ (p k ). Using the solution we derived in (2.35), one can explicitly check that Finally, in the limit β = 1, the Jack polynomials degenerate to Schur polynomials (characters of the un-deformed model) whose averages satisfy the analogous relation While we are not aware of an analytical proof of these relations, we have been able to check that they hold for all partitions of degree 9 and lower.

p = 2
The case of p = 2 is a generalization of the familiar Hermitean Gaussian matrix model, with generating function given by As we here require the n = −1 constraint in order to solve the model, we let ν = 0 and thus we can re-sum all the constraints starting from n = −1. In order to obtain an operator U whose diagonal is proportional to the dilatation operator D, we need to shift the weight of the re-summation as is an operator of degree 2, while L −1 is defined as in (2.27). As in the previous case D is the dilatation operator and it is of degree 0. An argument completely analogous to the one for p = 1 leads to the conclusion that U is a triangular operator with a 1-dimensional kernel, and therefore that the solution to the equation UZ p=2 (τ ) = 0 is unique up to normalization.
In order to give a W -algebra representation of the generating function we first consider the simpler case of a 1 = 0, then we re-introduce the parameter a 1 by shifting τ 1 ↦ τ 1 − a 1 . For a 1 = 0 we have the Gaussian case originally solved in [14] for which one can write Then we use the fact that together with the Virasoro constraint for n = −1 and a 1 = 0, to write the full solution as An explicit solution up to degree 3 is given by the correlators (2.46)

Averages of characters
As observed in [16], this model also satisfies the super-integrability property of characters.
Using the solution derived in the previous section one can check that the following relation holds Similarly, in the Schur limit where β = 1 we have which is also consistent with the result of [19] (as given in (2.19)) when a 1 = 0. These relations have been checked for all partitions up to degree 9.

Comments on p ≥ 3
Now consider higher values of p in (2.10). By re-summing all constraints in (2.26) for n ≥ −1 with weight (n + p)τ n+p we obtain the equation where the r.h.s. is a sum of shifted cut-and-join operators (2.50) and with deg W −p = p and deg K −k = k. The l.h.s. of (2.49) is of degree zero, therefore it corresponds to the diagonal part of the triangular operator U = ∑ n≥−1 (n + p)τ n+p U n , while W −p and K −k are of positive degree and therefore they represents the off-diagonal part of U.
For p ≥ 3 however, we immediately notice that the kernel of U is of dimension greater than 1. The kernel of the operator D − ∑ p−2 k=1 kτ k ∂ k is in fact infinite dimensional and corresponds to the span of all monomials which do not contain times τ k for k > p−2. For example, if p = 3 all monomials of the form τ 1 for all positive integer powers are annihilated by the diagonal part of U. This means that equation (2.49) does not provide a recursion relation expressing the corresponding correlator c {1,...,1} as a linear combination of correlators of lower degree. Consequently one should consider these coefficients as additional background data that needs to be specified independently in order to fully determine the generating function. As remarked in [25], for finite values of N one can always find additional relations between such correlators because only at most N of those can be linearly independent for a matrix of finite size. Therefore one can reduce the indeterminacy of the system of equations coming from the Virasoro constraints to a finite amount of information. Nevertheless, one cannot write a full solution for the generating function either in terms of correlators or in W -algebra representation.
We now present a formal way to repackage all the information that can be obtained from the recursion (for earlier attempts see [27][28][29]). From the integral representation of the generating function we can derive the additional identities where we have also explicitly written the dependence of the generating function on the coupling constants a k . If we substitute (2.52) in the l.h.s. of (2.49) we can rewrite the term − ∑ which is now no longer of zero degree in the times (in fact, since ∂ ∂a k has degree zero, every term in the sum has the same degree as τ k ) which means that it is not a diagonal operator. If we write W for the off-diagonal part of U, This constraint is still triangular (with respect to a basis of C[[τ 1 , τ 2 , . . . ]]) but now its diagonal component is the operator D, which we know has 1-dimensional kernel and in particular it is invertible over the complement of its kernel. A formal solution can now be obtained by splitting the generating function as and observing that D and (D − W ) are invertible operators when restricted to the image of W (which is contained in the complement of ker D), we obtain (2.58) This formal expression for the generating function automatically implements the additional constraints (2.52) but only for 1 ≤ k ≤ p − 2, precisely because in the operator W the derivatives with respect to a p−1 and a p do not appear. This implies that our solution (2.58) in general does not satisfy (2.52) if k = p − 1 or k = p. Without loss of generality then we can assume a p = 1 and a p−1 = 0 and repeat the argument that leads to (2.58). The final answer now is totally unambiguous and only depends on an appropriate choice of the correlation function (2.59) Because W contains derivatives in the variables a k , the recursion relations are no longer polynomial in the correlators c ρ (a). In fact we have that c ∅ (a) acts as a generating function for all the correlators that the recursion could not fix and the formal solution (2.58) expresses them as a-derivatives of c ∅ (a).
For example, if p = 3 we can write the solution up to degree 4 as .
It is curious to notice that while the integral representation of (2.59) is natural from the point of view of the definition of the matrix model, the solution of the Virasoro constraints does make sense also for an arbitrary function c ∅ (a 1 , . . . , a p−2 ) which does not necessarily admit an integral representation of that form.

Quantum models
We now shift our attention to the q-deformation of the classical models presented in the previous section. These will correspond to families of deformations depending on 1 or more parameters (typically q, t and in some cases r) which in the limit of those parameters going to 1 reduce the familiar examples already discussed. We refer to the q-deformed models as the quantum version of the Hermitean matrix model and to the degeneration limit q → 1 as their semi-classical approximation.
One more motivation for the name "quantum" is that we are able to identify such matrix models as the localized partition functions of certain supersymmetric quantum field theories in 3 dimensions. More explicitly, these correspond to theories with 4 supercharges placed on backgrounds of the form D 2 × q S 1 or S 3 b . As explained below, in some specific sense one can regard the D 2 × q S 1 partition function as a half of the partition function on S 3 b .
The q-deformation of the classical Virasoro constraints is a system of finite difference equations obtained by acting with operators satisfying a q-analogue of the Virasoro algebra. In the following, we mimic the derivation and solution of the Virasoro constraints in the q-case while simultaneously providing a detailed matching of the parameters between the quantum and the semi-classical case. Schematically, we have the identifications of parameters as in Table 1. Table 1. Matching of parameters between quantum and classical models. While the integer parameters N f and p can be straightforwardly identified, for the other parameters the identification is slightly less obvious. The β-deformation is obtained by identifying the adjoint mass as t = q β while ν can be related to the effective FI parameter through the parameter r. Similarly, the coupling constants a k are given by non-trivial functions of the masses u k . Quantum model Classical model

Definitions
We will now provide the details of the q-models which we intend to study, namely certain supersymmetric gauge theories in three dimensions.
We consider the partition function of a 3d N = 2 supersymmetric theory with gauge group U (N ) on the 3-manifold D 2 × q S 1 . More precisely, the geometry is that of a D 2 fibration over S 1 such that the D 2 fiber is rotated of a parameter q when going around the base. The U (1) holonomy q is identified with the quantum deformation parameter of the resulting matrix model. The partition function of an N = 2 theory on this geometry is sometimes referred to as the half-index [30] or holomorphic block [8].
Besides the N = 2 vector multiplet, we consider also an adjoint chiral multiplet of mass t and N f fundamental anti-chiral fields of masses u k . Moreover, we also turn on a Fayet-Iliopoulos (FI) parameter κ 1 ∈ C. The partition function can be computed via supersymmetric localization [8,9] with the result where the contour C is a middle dimensional cycle in (C × ) N defined by taking the product of N copies of the unit circle. Observe that there might be analytical issues with this naive choice of contour when the parameters are non-generic (see [8]). However we will avoid discussing such difficulties and assume that an appropriate contour exists by defining the generating function as an analytically well-defined solution to a set of partial differential equations obtained via algebraic manipulations of the integrand. More specifically, these equations will be the q-Virasoro constraints discussed in Section 3.2. In the generic case the two approaches are equivalent.
The integrand of the partition function is defined as the product of the classical contribution and the product of 1-loop determinants coming from the contributions of the vector, adjoint chiral and fundamental anti-chiral multiplets. Here (x; q) ∞ is the q-Pochhammer symbol defined in (A.1).
The partition function in (3.1) can then be interpreted as the matrix model where the measure can be seen as the q-deformed Vandermonde determinant, while the remaining contributions can be interpreted as q-deformations of the classical potential V (λ i ). The contribution Z cl D 2 ×qS 1 (λ) for instance has the form of a determinant insertion, however the actual exponential of the determinant will receive corrections from the measure ∆ q,t (λ). In order to make the connection to the potential in (2.16) we observe that the 1-loop determinant of the fundamental multiplets can be formally rewritten using the identity (A.4) as where p s (u) are the power-sum variables for the masses u k , Once we define the generating function by introducing the standard coupling to times {τ s }, we can see that (3.5) is a "potential" of the form which can be obtained by the shift of times .
Notice however that here all of the times must be shifted, as opposed to the classical case where only a finite number (corresponding to the integer p) had a non-trivial shift.
For a generic (polynomial) operator O = O(λ) we define its expectation value using the notation 9) where ⟨O⟩ N f is obtained by setting all the times to zero in the previous formula. Moreover, we assume that Z N f D 2 ×qS 1 (τ ) admits a formal power series expansion in times, whose coefficients are the correlators c ρ of the theory.
We pause here to explain the physical meaning of the generating function Z where ρ is the integer partition labeling the highest weight of the representation (see Appendix B for a short review of symmetric functions and Schur polynomials).

S 3 b
A different but intimately related model is that of an N = 2 YM-CS theory on the squashed 3-sphere S 3 b . We consider U (N ) gauge group and the same matter content as before. The 3d geometry is that defined by the equation where ω 1 , ω 2 ∈ R are the squashing parameters. The dependence of the partition function on the squashing is often indicated via a real parameter b such that b 2 = ω 2 ω 1 [10, [31][32][33]. We remark that, while geometrically it is natural to take the squashing parameters to be real valued, most of the formulas that we write in this paper are well-defined for arbitrary complex values. From now on, unless explicitly specified, we will assume ω 1 , ω 2 ∈ C.
Topologically, we can think of the 3-sphere as the gluing of two solid tori, i.e. two copies of D 2 × S 1 whose boundaries are identified via a modular transformation which acts by exchanging the two fundamental cycles of T 2 . Each half of the sphere can then be though to define a copy of the theory on D 2 × qα S 1 where now each copy has its own modular parameter q α , α = 1, 2, which we can express through the squashing parameters of the sphere as This simple geometric picture eventually leads to the very non-trivial property of factorization of the S 3 b partition function into a product of holomorphic blocks [34]. Here we are interested in yet another consequence of this factorization, namely the fact that the 3sphere partition function satisfies two independent sets of q-Virasoro constraints as shown in [13], hence the name modular double.
As in the previous model, the theory on S 3 b has an N = 2 vector, an adjoint chiral of mass M a and N f anti-chiral fundamental fields of masses m k . Again, we allow for a nonzero FI parameter κ 1 however in this case we are also forced to introduce a non-vanishing (bare) Chern-Simons (CS) level κ 2 . The reason for this is not physical in nature but rather it arises as a technical requirement necessary for having q-Virasoro constraints which can be written as PDEs in the time variables. In the case N f = 2 this was first shown in [17] where a unit CS level had to be introduced. Here we generalize that condition to arbitrary number of flavors N f ≥ 1 by imposing that 13) or equivalently, that the effective 5 CS level κ eff 2 ∶= κ 2 −N f 2 be vanishing. Observe that this condition is compatible with the cancellation of all perturbative anomalies even when the 5 In the presence of matter fields, the CS level receives quantum corrections, so that one can define an effective CS level κ eff 2 ∈ Z, which then has to satisfy a quantization condition for the theory to be free of anomalies. In particular, this implies that the bare CS level κ2 can be taken to be an half-integer number as long as we have an appropriate number of matter fields.
bare CS level is half-integral (i.e. N f is odd). We remark also that a non-zero effective CS level would correspond, from the point of view of the classical matrix model, to a potential term of the form − log 2 (λ), which would similarly spoil the derivation of the usual Virasoro constraints.
The partition function can be computed by means of supersymmetric localization techniques [10,11] and the result is given by where X i ∈ i R are Coulomb branch variables, Z cl which is written in terms of the double sine function S 2 (z ω) defined in (A.5). Analytical issues related to the convergence of the integral and its dependence on the physical parameters are addressed in Appendix D.
In order to simplify the notation we introduce the following set of exponentiated variables with ω = ω 1 + ω 2 . Introducing the complex parameter β = M a ω, we can also write t α = q β α .
The partition function in (3.14) then defines a quantum deformation of an eigenvalue matrix model with measure (3.18) and potential Observe that, for ω 1,2 in generic positions, one can use the identity (A.7) and regard the double sine function as the product of two q-Pochhammer symbols with arguments λ i,1 and λ i,2 , respectively. This can be used to argue that the S 3 b generating function satisfies two independent set of q-Virasoro constraints, however it can be shown (see [17]) that this latter property holds for any value of the squashing parameters, even when (A.7) does not apply.
Because of the presence of two separate (but not independent) sets of variables {λ i,α }, we can define a "doubled" generating function by coupling each set to its own copy of auxiliary time variables, where now {τ s,α } act as conjugate variables for the power-sum observables As a formal power series in times, the generating function can be written as for ρ, ρ ′ integer partitions of arbitrary sizes and are the correlation functions of products of the power-sums.
Expectation values of Wilson loops are defined equivalently to the case of D 2 × q S 1 given in (3.10), however they now carry a dependence on α due to the power-sum variables (3.21) in the argument of the Schur polynomial. Physically, this corresponds to the fact that on S 3 b there are two BPS cycles, one for each solid torus, and therefore the holonomy of the gauge connection must carry an index α distinguishing between the two.
It is now worth noting that we can formally obtain the generating function on D 2 × q S 1 (3.7) from the generating function on S 3 b (3.20) by setting one set of times to zero, {τ s,2 = 0} for instance, (3.24) These generating functions are then equivalent in the sense that they satisfy the same set of constraints as will be shown below in (3.34) and (3.41), however they do not have the same integral representation. Thus, whenever we encounter ambiguities in defining the contour C in the partition function of D 2 × q S 1 in (3.1) we can resolve them by choosing any contour which is consistent with (3.24). In other words, we can choose any contour such that the coefficients in the power series expansion of Z

q-Virasoro constraints
Mirroring the procedure of the classical case, we will now derive constraint equations for the gauge theory on D 2 × q S 1 and then briefly summarize the corresponding procedure in the case of S 3 b . The techniques used to derive the constraints are those developed in [12,17] therefore we avoid presenting each step of the computation. The most important differences with respect to those cases are the presence of an arbitrary number of flavors N f and the introduction of the additional deformation parameter r.
We now show that the generating function (3.7) satisfies a set of first order q-difference equations which take the name of q-Virasoro constraints. In order to derive these q-Virasoro constraints we introduce the finite difference operatorM i defined aŝ for f a function of the gauge variables {λ i }. The constraints are obtained by substituting in (2.23) the partial derivative ∂ ∂λ i with the difference operator 6 The vanishing of the q-variation of the generating function can then be written as with ". . . " denoting the integrand.
For the sake of simplicity of the computation, we sum all such equations for the integer n ranging over all of Z and multiply each component by the corresponding power z n of an auxiliary formal variable z. We finally end up with the single q-Virasoro constraint equation and each independent constraint can be recovered by expanding in the given power of z. 6 Observe that in the semi-classical limit t → 1 and q → 1, the function Gi(λ; t) goes to 1 while The LHS and RHS can be computed independently. Starting with the RHS we have where we recall the definition (3.6) of the power-sum variable p s (u), while (3.30) For order n ≥ 0 in the variable z all the insertions contain only positive powers of λ i , therefore the corresponding constraints can be expressed as PDEs in the time variables. At order n = −1 in z, there are negative power contributions in λ i which must vanish for the constraint to be meaningful. The offending term is which vanishes precisely if we let t = q β with β ∈ C and impose the balancing condition which is consistent with the choice of FI parameter in [13]. More generally we define and we refer to ν as the balancing parameter, which vanishes exactly when the balancing condition is satisfied. The q-Virasoro constraints can then be written as (3.34) The precise relation between these constraints and the generators of the q-Virasoro algebra are delucidated in Appendix C.
For the purpose of actually writing these constraints as recursion relations for the correlation functions c ρ , it is more convenient to rewrite the shift of times as multiplication by a polynomial in z −1 and the masses u k , using the identity where the coefficients A k = A k (u) are defined as antisymmetric Schur polynomials in the masses u k , Notice that antisymmetric Schur polynomials of degree higher than the number of variables {u k } are identically zero, therefore we have only a finite number of coefficients A k and the series expansion in (3.35) is truncated to degree N f .

S 3 b
Let us now derive the constraints for the S 3 b generating function (3.20). As already mentioned, there are two sets of independent q-Virasoro constraints that we can write. They correspond to the Ward identities obtained by acting with q-difference operators that shift each set of variables {λ i,α } separately. Namely, we can define the operatorsM i,α which act as:M (3.37) Because of the periodicity of the exponential function, it follows that on the exponentiated variables {λ i,α } the shift acts multiplicatively aŝ which explains why we use the same symbol as in (3.25). Proceeding as above, we define functions G i,α (λ; t α ) as and compute the Ward identities as the integral equations Keeping α fixed, we sum over all n ∈ Z with weight z n and obtain the q-Virasoro constraint (3.41) where we defined the balancing parameters ν and r α as If we explicitly expand (3.41) in powers of z, we have that for n ≥ 0 all equations are free of expectation values of negative powers of {λ i,α } and therefore can be rewritten as PDEs for the generating function. At order n = −1, there are negative power contributions which must be cancelled. An explicit computation shows that these offending terms vanish exactly when ν = 0 (so that r α = 1). We call this the balancing condition. Observe that ν corresponds to the combination √ βα of [13].
Finally, by applying the reasoning of [17], we find that the correlators c ρ;ρ ′ in (3.23) obtained as solutions to the above constraints, do factorize according to where c ρ;∅ is a correlation function containing only the {λ i,1 }, and c ∅;ρ ′ is a correlation function only containing the {λ i,2 }. The factorization of the correlators then implies that (3.22) can be written as a product where each half is a formal power series in one set of times only, and it satisfies its own copy of the q-Virasoro constraints. Here we also used that the empty correlator is the partition function, i.e. the generating function evaluated on all the times equal to zero, (0). Using this fact, we deduce that up to normalization, the correlators of the S 3 b theory are products of correlators of D 2 × q S 1 theories. The formula (3.44) gives a precise meaning to the equivalence proposed in (3.24).
Remark. As observed in the pole analysis of the partition function in [17,Appendix C], the choice of shift operatorM i,α is motivated by the requirement that when we shift the contour of integration in the LHS of the constraint equation in order to reabsorb the action of the shift in (3.37), we should not cross any poles of the integrand. Since the only poles that can come in between the two contours are the ones due to the fundamental matter fields, the choice of the chirality of these fields is restricted by the sign of the shift in the variables. Namely, for negative shift of the X i variables, we need to shift the contour to the left and therefore we cannot have any poles to the left of the imaginary axis. This implies that all fundamental matter fields should be anti-chiral rather than chiral. Their precise number is then fixed by the value of the bare CS level as in (3.13). However, if we were to use the opposite shift operatorM −1 i,α , we would need to shift the contour in the opposite direction, and by doing so we might cross some poles along the way. To avoid this problem we can substitute all anti-chiral fundamental multiplets with chiral fundamentals. The change of chirality inverts the position of the poles of the corresponding double sine functions thus allowing us to shift the contour. By reproducing the analysis of the q-Virasoro constraints we find that the corresponding condition on the number of fundamental chiral fields is N f = −2κ 2 , which can only be satisfied for negative bare CS level. The rest of the derivation is completely equivalent provided we adopt the substitution of parameters This is indeed compatible with the symmetry properties of the partition function proved in [35,Proposition 5.3.16].

Solution of the constraints
We now study the conditions under which the q-Virasoro constraint equations admit solutions and whether these solutions are unique or not. As we shall show, the only cases in which the solution is uniquely defined are those of N f = 1, 2 and in such cases we use the solution (computed algorithmically up to a certain finite order) to explicitly verify the property of averages of Macdonald polynomials conjectured in [19]. Having found an explicit solution, we study the semi-classical limit and propose an identification with the models of Section 2.
Before actually solving the constraints (3.34) and (3.41), we observe that because of the factorization property (3.44), we just need to find a solution to the q-Virasoro constraints on D 2 × q S 1 . Up to normalization then, a solution for the generating function on S 3 b can be obtained by taking the product of correlation functions as in (3.43). For this reason in this section we consider an abstract generating function Z N f (τ ) which depends on a single set of times and satisfies one set of q-Virasoro constraints with respect to those times. The first step in finding a solution to (3.34) is to rewrite the whole set of constraints as an equation where U q,t is the operator obtained by appropriately re-summing the individual constraints, similarly to how we did in the classical case. What we then find is that U q,t is a quantum deformation of the operators in (2.29) and (2.40). Moreover, we observe that the q-deformation preserves the triangular form 7 of the operator. In order to see this, we manipulate (3.34) by separating the exponential of the times from the exponential of the p s (u) and then applying (3.35) to obtain (3.47) The operator U q,t is defined by expanding in powers of z and re-summing over all n ≥ −1 with weight (n + N f )τ n+N f . Finally, one can check by explicitly writing U q,t as a linear operator on the formal ring C[[τ 1 , τ 2 , . . . ]], that it assumes the form of a semi-infinite triangular matrix 8 .
If we split the constraint operator as where D q,t is the diagonal part and W q,t is the strictly upper triangular, then as we will show in the following sections, we can write up to some overall coefficient which does not depend on the times. For N f = 1, 2 the operator is diagonal with all eigenvalues different from zero except for the one associated to the constant monomial (labeled by the empty partition), hence we deduce that both D q,t and U q,t have a 1-dimensional kernel just like in the classical case of p = 1, 2. This means that the solution to the equation (3.46) is unique up to normalization of the empty correlator c ∅ . For higher N f the kernel becomes infinite dimensional and the solution is no longer uniquely defined. This suggest a strong parallel with the classical models upon identification of the parameters N f and p. In the following sections we analyze this question in detail and provide a concrete limiting procedure which explicitly shows the correspondence. 7 The main algebraic difference between Uq,t and U classical is that the former is only triangular with respect to the full basis of monomials with the total ordering induced by integer partitions, while the latter is also triangular with respect to the partial order given by the monomial degree (size of the partition). 8 Observe that every triangular system of linear equations is equivalent to a set of recursion relations between the components of the solution. In our case it corresponds to an infinite set of finite-step recursion relations between the correlators cρ.

N f = 1
We first consider the case of N f = 1. Using Cauchy's formula (B.8) we can expand (3.47) in powers z n for n ≥ 0 to get (3.50) If we define the differential operator U n as the operator such that U n Z N f =1 (τ ) = 0 is the n-th q-Virasoro constraint (3.50), then we re-sum all such operators over n ≥ 0 to define which we can write explicitly as where the first term in the right hand side is an operator of degree 0 while the other two are of degree 1. The degree zero part is the one that contains the diagonal operator D q,t which is defined by expanding the symmetric Schur as in (B.4) and picking only the term p n n (all other terms are higher derivatives in the times and therefore cannot be diagonal), The remaining terms in the r.h.s. of (3.52) give the definition of −W q,t . Now we have that the full set of q-Virasoro constraints is equivalent to the statement that Z N f =1 (τ ) is in the kernel of the operator U q,t . The fact that this operator is not homogeneous in degree then corresponds to a set of recursion relations between correlators in degree d and those in degree d − 1. In particular, the recursion will allow us to write c ρ as a linear combination of all c ρ ′ such that deg(ρ ′ ) ≥ deg(ρ) − 1 and ρ ′ < ρ. Because the linear operator in (3.52) has a 1-dimensional kernel, the generating function Z N f =1 (τ ) is unique up to a choice of a normalization constant c ∅ ≡ Z N f =1 (0). An exact solution up to degree 3 is then given by the correlators 54) which are rational functions in the parameters q, t, r and A 1 . Moreover, we have that all correlators of degree higher than 1 are proportional to c {1} precisely as in (2.36).
Observe that the dependence on the rank of the gauge group N comes only through the powers of t, and in the large N limit the correlators behave as where ρ is the size of the corresponding partition.

Semi-classical limit
As we already observed, the case of N f = 1 bears many similarities with the classical model with p = 1. The precise relation can be established by a semi-classical limit procedure in which we introduce a perturbation parameter ̵ h such that q = e ̵ h . In the limit ̵ h → 1 both the constraint equations and their solution match exactly with the formulas for the classical case p = 1. The matching goes as follows.
First we consider the measure ∆ q,t in the generating function. Assuming t = q β with β ∈ Z, we can write which in the limit ̵ h → 0 (or q → 1) clearly gives a β-deformed Vandermonde determinant ∆ 2β (λ) multiplied by a determinant insertion of power −β(N − 1). Combining this determinant with the one coming from the FI term and the one in the measure d λ i λ i , we have where we used (3.33) to write the exponent as the balancing parameter ν.
The last terms to match are the 1-loop determinant of the fundamental anti-chiral and the potential V (λ i ). From (3.8) we see that the shift of times is singular in the limit q → 1. In order to make this limit well defined, we also scale the mass u 1 as with a 1 a positive constant independent of ̵ h. This way we can re-write the 1-loop determinant as Hence we obtain the same integrand as in (2.14). Notice that in order to get the correct Vandermonde term in the semi-classical limit, we had to assume that β is an integer number. For arbitrary complex values however (3.56) does not make sense, nevertheless we know that the generating function admits an analytic continuation in β, then we can assume that our limit still makes sense even for non-integer complex numbers.
With the choice of mass as in (3.58) the semi-classical limit q → 1 is also well-defined at the level of the q-Virasoro constraints. To this end we assume the following power series expansion in the parameter ̵ h around 0 both for the generating function and for the constraint operators themselves, where the first non-trivial equation comes by setting the coefficient of ̵ h to zero. We thus have that Z This is precisely the operator in (2.28). Hence we have the identification of generating functions We remark that in the case of the generating function on S 3 b there are two sets of exponentiated masses u 1,α but only one actual fundamental mass m 1 . Therefore our parametrization in (3.58) leads to a well-defined limit only for one copy of q-Virasoro at a given time. Namely, if the mass has been chosen so that one copy of the constraints becomes usual Virasoro in the semi-classical limit, then the other copy is not well-behaved in that limit.

Averages of Macdonald Polynomials
From the point of view of the q-deformed matrix models it is natural to consider expectation values of Macdonald polynomials. These polynomials are the natural quantum generalization of the Schur polynomials and as such they can be interpreted as the correct quantum characters of the model. We now show that their averages do indeed satisfy a special property of the form ⟨character⟩ = character.

Introducing the functions
we have the identity We have checked that this formula holds for all partitions up to degree 6. In principle our solution should allow to check up to arbitrary finite order, however for higher degrees the computation becomes quickly too impractical even for computer calculations.

N f = 2
In order to show that a solution in this case exists and is unique we need to use all constraints for n ≥ −1 and in particular we need the n = −1 constraint to be well-defined. This implies the balancing condition r = 1, (i.e. ν = 0). The n-th q-Virasoro constraint operator U n can be written as so that the quantum operator U q,t = ∑ ∞ n=−1 (n + 2)τ n+2 U n takes the form where the operator in the first line is of degree 0, the one in the second line is of degree 1 and those in the other three lines are of degree 2.
Similarly to the previous case, the operator U q,t is triangular and it has a 1-dimensional kernel, spanned by the solution Z N f =2 (τ ) of the recursion. Again, U q,t can be split as (3.48) into a diagonal part and an off-diagonal part which we collectively call W q,t . While D q,t has degree zero, W q,t is not homogeneous and it contains terms of degree 0, 1 and 2. The recursion in this case has a longer but still finite step and an exact solution can be computed from the initial condition c ∅ by solving the equation order by order in the times {τ s }. For concreteness we present the solution in terms of correlators up to degree 3, 74) which is consistent with the result in [17] (if we identify A ≡ A 1 and B ≡ A 2 ). We observe that the correlators are rational functions of the parameters q, t, A 1,2 , and that in the large N limit they behave as (3.75)

Semi-classical limit
In order to make the semi-classical limit q = e ̵ h with ̵ h → 0 well-defined at the level of the shift of times in (3.8), we need to appropriately choose the masses u 1 , u 2 so that they scale non-trivially with ̵ h in such a way that the shift is finite. By imposing the conditions for a 1 , a 2 two (positive) constants, we can parametrize the choice of masses as (3.77) up to permutation u 1 ↔ u 2 . The anti-chiral 1-loop determinant of the fundamental antichirals becomes where the coefficients of the higher powers of λ i are completely determined by a 1,2 and, more importantly, they all vanish in the limit q → 1. With this parametrization we also have that the variables A 1 , A 2 behave as Expanding the operator U n in powers of ̵ h as in (3.62), the first non-trivial contribution is which we immediately recognize as the Virasoro constraint operator for the classical matrix model in (2.39). By re-summing with weight (n + 2)τ n+2 over n ≥ −1, we get the relation − U

Averages of Macdonald Polynomials
Evaluating the average of Macdonald polynomials on the explicit solution that we found, we are able to write the following identity Assuming that this formula holds for all partitions ρ, we can compute the semi-classical limit t = q β and q → 1 with the choice of masses in (3.77). The result matches exactly with formula (2.47) for the average of Jack polynomials at p = 2.

Comments on N f ≥ 3
Our analysis indicates that the only cases that admit a full and unique solution (independent of the normalization of the empty correlator) are those of N f = 1, 2. For higher values of N f we have a similar situation to that of the classical models at p ≥ 3. The q-Virasoro constraints have an infinite dimensional kernel, corresponding to the fact that one needs to specify more initial conditions to solve the recursion. We believe that the situation here mirrors what we discussed in Section 2.3.3, however the actual formulas are too cumbersome and unilluminating to write down explicitly. We remark however that computer calculations suggest that the actual recursion relations can be solved iteratively by the same procedure delineated there. We do however expect that a more subtle study of possible analytical issues is required in this q-deformed case, especially with regards to the dependence on the mass parameters u k . The definition of the matrix model for instance might have ambiguities related to the choice of contour, similar to the so called Dijkgraaf-Vafa phase of a classical matrix model [36,37].
The semi-classical behavior is also not so straightforward. This limit can, in fact, be defined by the choice of masses which solves the equations In terms of the u k , these form a system of N f equations of increasing degree up to N f . Over the complex numbers there are N f ! solutions which are all equivalent upon permutations of the u k , however writing such a solution explicitly is in general not possible. What we expect is that upon substitution in (3.8) we get a shift of the times such that only the first N f terms are non-vanishing in the limit q → 1. Then the constants a k parametrizing the solution, can be identified with the coupling constants of the semi-classical model where p = N f .

Conclusion
In this paper we give an outline of the procedure to obtain and recursively solve the Virasoro and q-Virasoro constraints in the case of the classical β-deformed Hermitean 1-matrix model and the quantum 3d N = 2 theory with U (N ) gauge group on D 2 × q S 1 and S 3 b . We present the solution of the models in terms of explicit expressions for the first few correlators, and additionally in the classical case we can also express the solution using the W -representation of the generating function. Moreover, we deduce novel formulas of the form ⟨character⟩ = character in the spririt of [19] as given in (2.37), (2.47), (3.69) and (3.82). Finally, we explicitly match the classical models with their q-deformations by showing the existence of a well-defined semi-classical expansion around a formal deformation parameter q = e ̵ h , for small ̵ h.
There are several directions for further investigation.
• One obvious direction is to extend the W -representation of the generating function to the quantum models, an investigation which has been initiated in [26]. This is a much more difficult task than in the classical case due to the more involved constraint equations.
• Another interesting direction might be to try to find an analytical proof of the formulas for expectation values of characters. These were all unexpected a priori, and the formulas appearing in this paper have only been verified by explicitly checking the equations for partitions up to some finite degree.
• Another open question is regarding the condition (3.13) on the effective CS level. In particular how one can treat models when κ eff 2 is different from zero.
• Furthermore, another question is whether the procedure of obtaining and solving the Virasoro constraints can be generalized to other root systems, in the sense that the above derivation is valid only for a U (N ) gauge group or A N from a matrix model perspective. The question is then if there are other versions of the derivation for other Lie groups, which then both offers a definition of the matrix model, together with a set of Virasoro constraints which may or may not be solvable.
• Another interesting question to consider, is why the system of equations in the Virasoro constraints are triangular and give rise to a finite recursion relation. From the explicit Virasoro constraints it is clear that the equations are triangular, but why this is so in the first place is not obvious.
• There is also the question if there is a deeper reason to why p ≥ 3 in the classical case and N f ≥ 3 in the quantum case cannot be solved. At the level of the constraint equations we understand why this is not solvable, but if there is a more fundamental reason to why this is the case is still not clear.
When z < 1 and q ≠ 1 one can define the quantum dilogarithm [38] which can be used to rewrite the q-Pochhammer as Secondly, we introduce the double sine function [39,40]. For ω ≡ (ω 1 , ω 2 ) ∈ C 2 with Re(ω 1 ) > 0, Re(ω 2 ) > 0 and z ∈ C, the double sine function is defined by the regularized infinite product which satisfies the inversion property For generic values of the parameters ω 1 , ω 2 such that Im( ω 2 ω 1 ) ≠ 0, the double Sine function has the following infinite product representation where B 22 (z ω) is the double Bernoulli polynomial

B Symmetric functions and characters
The Schur polynomials denoted by Schur γ (λ 1 , . . . , λ N ) are labeled by partitions γ and are defined as the irreducible characters of the group U (N ). As such they form an orthonormal (linear) basis in the space of all polynomial characters which, by definition, are invariant under the action of the Weyl group S N . This means that Schur polynomials form a basis for all symmetric functions.
For any partition γ = {γ 1 , . . . , γ N } whose elements obey γ 1 ≥ ⋅ ⋅ ⋅ ≥ γ N ≥ 0 one can compute the Schur polynomial Schur γ (λ 1 , . . . , λ N ) as a ratio of determinants as in the Weyl character formula Introducing the power-sum variables p k = ∑ N i=1 λ k i one can expand Schur polynomials as where χ γ ρ is the character of the representation of the symmetric group indexed by the partition γ evaluated at elements of cycle type ρ. Here we also introduced the notation Aut(ρ) for the automorphism group of the partition ρ, namely the group of permutations of the parts of ρ which are of equal length. The order of this group can be computed as For symmetric Schur polynomials (B.2) takes the simpler form where γ ⊢ m denotes that γ is an integer partition of m.
If we define the scalar product of two symmetric functions by The Jack and Macdonald polynomials are defined as 1-or 2-parameter deformations of the Schur polynomials (see [41]). More specifically, for any β ∈ R >0 we define Jack polynomials Jack γ (p k ) as the symmetric functions orthogonal with respect to the inner product Given parameters q < 1 and t < 1 we can define Macdonald polynomials Macdonald γ (p k ) as a family of symmetric functions orthogonal with respect to the inner product For concreteness we only consider Jack and Macdonald polynomials in the P -basis. In the limit t → q, Macdonald polynomials degenerate to Schur polynomials, while for t = q β they give the Jack polynomials.

C Relating constraint generators to q-Virasoro generators
We now wish to obtain a relation between the generators of the constraint U n in the quantum case and the generators of the q-Virasoro algebraT n . Similar to the analysis in [17], one can introduce the function ψ(z) defined as ψ(z) = q −1 2 t 1 2 r −1 2 exp ∞ s=1 z −s (1 − q s ) (1 + q s t −s ) τ s (C.1) and the generator currentT (z) = ∑ n∈ZTn z n aŝ Here the only differences compared to the standard generator current T (z) of the q-Virasoro algebra [42], are therefore the factors of r ±1 2 and the shift of times proportional to p s (u).
We can also denote the eigenvalues of the generatorT n for n ≥ −1 aŝ That the operatorT n is diagonal also for n = −1 is not obvious, although we show below that this is the case.
Then, we can rewrite the left hand side of the constraints in (3.34) as so that one can read off the eigenvalue for theT 0 by equating (C.5) and the right hand side of (3.34) to find For theT −1 operator, the eigenvalue is only well defined in the case of ν → 0 or r = q ν → 1 (when we also recover equation (3.37) of [17]), in which case it iŝ From the dependence on p 1 (u) in the above, we can also see that the eigenvalue of T −1 vanishes.
We now recall that the q-Virasoro constraints were written using the generators U n as U n≥0 Z N f D 2 ×qS 1 (τ ) = 0. Then the connection between the generators U s and the generator currentT (z) can be written as Let us look at the fundamental anti-chiral contribution to the integrand. These are of the form S 2 (−X i − m k ω) −1 = S 2 (X i + m k + ω ω) (D.2) where we used the inversion formula (A.6). Plugging z = X i + m k + ω in (D.1) we obtain the asymptotic approximation Therefore, at large X i the fundamental chirals behave as shifts in the CS level and FI parameter. Putting together all the potential terms we get e − π i κ 2 ω 1 ω 2
• For Im(X i ) > 0 (i.e. = +1), 2κ 2 − N f = 0 and the quadratic term X 2 i vanishes. The linear term becomes dominant and the exponential goes to zero at infinity if Re (ω(ν + 1)) = Re ⎛ ⎝ In our conventions ω ∈ R >0 , so that this constraint is equivalent to the requirement Re(ν) > −1 which corresponds to the condition for convergence of the integral that we observed in the classical model in (2.14). For N f > 1 we have ν = 0 so that (D.7) is trivially satisfied.
This analysis then implies that the integral can be computed by closing the contour in the right-hand-plane and taking the residue as in Figure 2.  Physically speaking, one can map the real part of the complex masses to the R-charge 9 In order to see this, one should write ∆ S (X) as the product ∏ α=1,2 ∆q,t(λα) and then apply (3.56).
of the corresponding field as follows: Re(m k ) = − ω 2 R k and Re(M a ) = ω 2 R a , with R k the R-charge of the k-th fundamental field and similarly R a is the R-charge of the adjoint field, while Re(κ 1 ) = ω 2 R monopole is the R-charge of the bare monopole operator [44]. Rewriting (D.7) in terms of the R-charges of the fermions and using the fact that the gauginos carry R-charge 1, we have so that convergence of the integral imposes the positivity of the effective R-charge of the monopole operator, while the balancing condition (ν = 0) further restricts this charge to be 2. This considerations suggest that there might be non-perturbative effects giving rise to a monopole superpotential when ν = 0 similar to those studied in [45].