Exact SUSY Wilson loops on $S^3$ from $q$-Virasoro constraints

Using the ideas from the BPS/CFT correspondence, we give an explicit recursive formula for computing supersymmetric Wilson loop averages in 3d $\mathcal{N}=2$ Yang-Mills-Chern-Simons $U(N)$ theory on the squashed sphere $S^3_b$ with one adjoint chiral and two antichiral fundamental multiplets, for specific values of Chern-Simons level $\kappa_2$ and Fayet-Illiopoulos parameter $\kappa_1$. For these values of $\kappa_1$ and $\kappa_2$ the north and south pole turn out to be completely independent, and therefore Wilson loop averages factorize into answers for the two constituent $D^2 \times S^1$ theories. In particular, our formula provides results for the theory on the round sphere when the squashing is removed.


Introduction
During the last 30 years there has been a vast development in our understanding of supersymmetric gauge theories in various dimensions. For many supersymmetric gauge theories the partition functions and the expectation values of certain protected BPS observables can be calculated exactly. It has been observed that these exact gauge theory quantities can be expressed in terms of two dimensional conformal field theories (or their deformations), a phenomenon known as BPS/CFT correspondence [1,2]. One famous example of the BPS/CFT correspondence is the AGT correspondence which relates the 4d S-class theories to 2d Liouville and Toda models [3,4]. By now we know many more concrete examples of the BPS/CFT correspondence. In this paper we concentrate on the concrete application of the BPS/CFT correspondence to 3d supersymmetric gauge theory. In particular, we will show how this correspondence leads to explicit formulas for the expectation values of supersymmetric Wilson loops.
Starting from the work [5], there has been many explicit calculations of the partition functions and other BPS observables for supersymmetric gauge theories on compact manifolds, see [6] for a review. In these calculations the main tool is equivariant localization in the space of fields, and the final answers are typically expressed in terms of finite dimensional matrix models which are generically rather complicated. In the case of the 3d N = 2 Yang-Mills-Chern-Simons (YM-CS) theories, the corresponding matrix models were derived in [7] for the round sphere S 3 and in [8,9] for the squashed sphere S 3 b . The expectation value of the supersymmetric Wilson loop corresponds to the specific insertion of a Schur polynomial in the matrix model and it is convenient to combine them into the generating function Z (τ 1 , τ 2 ). For a generic value of the squashing parameter b, this generating function encodes all Wilson loops in arbitrary representations. In [10] it has been observed that Z (τ 1 , τ 2 ) for 3d N = 2 YM-CS U (N ) theory coupled to adjoint and possibly (anti)fundamental chiral multiplets has a free field representation in terms of vertex operators and screening charges of the q-Virasoro modular double (this observation is based on earlier works [11] and [12,13]). Upon fixing some parameters, the generating function Z (τ 1 , τ 2 ) satisfies two commuting sets of q-Virasoro constraints which provide the Ward identities for the corresponding matrix model. However, this free field construction is formal and it breaks down in the case of the round sphere b = 1. In this paper we would like to address the numerous analytical issues and solve these Ward identities explicitly.
The present paper is the development of the ideas and techniques of [14], where the YM-CS living on D 2 × S 1 was worked out in detail. There the Ward identities (the q-Virasoro constraints) for the matrix model were derived by inserting certain q-difference operators under the integral with some analytical issues being addressed, and finally the Ward identities lead to the iterative solution for all correlators in the corresponding matrix model. Here we consider a particular supersymmetric gauge theory -the N = 2 U (N ) Yang-Mills-Chern-Simons (YM-CS) theory coupled to one adjoint and two fundamental anti-chiral matter multiplets on the squashed sphere background S 3 b (see Section 2 for the precise definition). Going to the squashed sphere case requires, among other things, a new careful analysis of the poles coming from contributions of gauge and matter multiplets (outlined in Appendix C). This analysis needs to be performed case by case for every theory one considers. How generic the class of theories for which our procedure works is therefore a subject of further research. At the end we provide a simple and efficient way to algorithmically calculate (supersymmetric) Wilson loop averages in any concrete representation. This procedure, which can be readily cast into computer program form, is the main new contribution of this paper. As a rough illustration of our result, for a Wilson loop around the north (α = 1) or south (α = 2) pole in representation R = {1, 1}, i.e. the rank 2 antisymmetric, we get , (1.1) where A, B, and t α are functions of the fundamental and adjoint masses and the squashing parameters (see Sections 2 and 3). In this paper we will always be discussing normalized expectation values of Wilson loops.
Furthermore, even for the particular theory on S 3 b we managed to get the whole scheme of [14] working only for special values of the Chern-Simons (CS) level κ 2 and Fayet-Illiopoulos (FI) parameter κ 1 (see Section 3 for details). While these restrictions on κ 1 and κ 2 appear to be purely technical, i.e. they are needed to drastically simplify parts of the computation, at the moment it is not clear how to lift them. Therefore, some conceptual underlying reason for these restrictions may exist. We hope to address and push these limitations in the future.
In addition to their practical 3d gauge theory use, the present analysis and the corresponding Ward identities (3.28) are interesting from a purely matrix model point of view as well. They are nothing but q-Virasoro constraints (see Section 3.3), where the choice of representation of the Heisenberg generators depends on the adjoint and antifundamental masses. This was anticipated already in [10], but in this paper we pay careful attention to various analytical issues, which required introduction of antifundamental multiplets. Thus, it puts the formal derivation of [10] on a firm footing. Interestingly, as one takes the round sphere limit (Section 5), the representation of the Heisenberg generators becomes singular and the free field representation fails. However, the Ward identities admit the well-defined limit and our result is still valid for round S 3 . It can be noted that the round sphere limit does not correspond to the standard semi-classical limit of the q-Virasoro algebra collapsing to (usual) Virasoro algebra.
There is yet another angle from which our work may be interesting. Namely, a certain q− (but not t) deformed matrix model (the BEM-model [15]) which calculates colored HOMFLY polynomials for torus knots. Generalizations of this model, both in the direction beyond torus knots and in the direction of turning on the t parameter (i.e. going to Khovanov, Khovanov-Rozansky and superpolynomials) are much sought for. Once such (q, t)-deformation would be available, it would be very interesting to see how the techniques developed here, help to explain certain strange phenomena of the (q, t)-world such as chamber structures [16] and nimble evolution [17].
The paper is organized as follows. We continue in Section 2 with the definition of the gauge theory that we consider. In Section 3 we derive the q-Virasoro constraints using the insertion of a certain finite difference operator, and we also interpret the result using the free field representation. We then recursively solve the constraints and obtain explicit expressions for the expectation values of the first few supersymmetric Wilson loops in Section 4. In Section 5 we take several interesting limits of the result. Finally we summarize and suggest directions for further research in Section 6. Details of Schur polynomials and partitions, special functions and the difference operator and the shift of integration contour are left to the appendices.

Gauge theory on the squashed 3-sphere
In this section we give the definition of the theory that we will be working with and we also review some technical aspects regarding partition functions of 3d N = 2 YM-CS gauge theories.
We are interested in theories with unitary gauge group U (N ) and (anti-)chiral matter in the fundamental or adjoint representations. Such theories can be placed on curved compact backgrounds while still preserving 2 supercharges as shown in [7][8][9] and [18][19][20][21][22]. We work on a compact manifold of the form of a squashed 3-sphere S 3 b , defined as the locus where ω 1,2 are the (real) parameters of the squashing. It will also be useful to define the combination Upon analytic continuation of the partition function we can take the squashing parameters to be arbitrary non-zero complex numbers.
Another useful way to describe the squashed sphere background is that of a singular elliptic fibration over an interval. In this picture, one of the two cycles of the torus fiber shrinks to a point at one edge of the base while the other cycle shrinks to a point at the opposite edge.
If we cut open the interval at its midpoint, then the restriction of the fibration over each of the two smaller segments has the topology of a solid torus D 2 × S 1 with the degenerate fiber identified with the locus {0} × S 1 (where {0} is the center of the 2-disk). The gluing along the boundary ∂(D 2 × S 1 ) ≅ S 1 × S 1 is done via a diffeomorphism that exchanges the two fundamental cycles of the torus. Using the intrinsic coordinates θ, φ, χ given by we can identify θ ∈ [0, π 2 ] with the coordinate along the base and φ, χ ∈ [0, 2π] with the coordinates on the torus fibers. The singular fibers are then given by the cycles at θ = 0 and θ = π 2 which are parametrized by φ and χ, respectively. Upon using supersymmetric localization one finds that the 1-loop contributions of the gauge and matter multiplets are given in terms of products of double sine functions 1 S 2 (z ω) as summarized in Table 1.
1 The comparison between our notation and the literature is that S2 (ω 2 − iX ω) = s b (X), with ω1 = Here α ∈ ∆ are the roots of the algebra (we always exclude the zero root) and w are the weights of the representation R, while M R are the masses of the (anti)chiral fields. With X = (X 1 , . . . , X N ) we indicate the gauge variables, i.e. the integration variables of the localized partition function, taking values in the Cartan of the gauge group. For the case of U (N ) the roots are differences of fundamental weights w i so that we can write where the X i are imaginary numbers.
In addition to the 1-loop determinants we also allow for a CS term with level κ 2 ∈ Z and, since the gauge group has a U (1) in the center, an FI term with complexified parameter κ 1 ∈ C.
For technical reasons which will become clear in the following sections, we further specialize to a theory with a single U (N ) gauge group together with 1 adjoint massive chiral and 2 fundamental antichirals with masses µ, ν ∈ C. The partition function can then be written explicitly as the integral 2 (2.7) 2 Up to an overall multiplicative factor, this partition function coincides with the "level 6" integral Notice that the choice of CS level there is compatible with the one we have in Section 3.
Following the mathematical literature we denote the combination of the vector's and adjoint chiral's 1-loop determinants as the function ∆ S (X) which from the point of view of the matrix model theory represents a generalization of the Vandermonde determinant of U (N ). For more details on this see [24] and references therein.
Of great importance for the gauge theory is the computation of expectation values of supersymmetric Wilson loop operators. These correspond to quantum averages of traces of the holonomy of the gauge connection around some supersymmetric closed curves inside the spacetime manifold. Such supersymmetry preserving loops are referred to as 1 2 -BPS loops and, for generic ω 1 , ω 2 , are exactly the two singular fibers at θ = 0 and θ = π 2 . Whenever the ratio of the squashing parameters is a rational number we also have a second family of 1 2 -BPS loops at θ ≠ 0, θ 2 and wrapping the regular fibers according to the equation ω 1 φ+ω 2 χ = const [25]. By definition these cycles are torus knots inside of S 3 b . In this paper we will only consider the case in which the Wilson loops wrap one or both of the singular fibers. Concretely, the traces are taken over arbitrary irreducible representations R ρ of the gauge Lie algebra which for U (N ) are given as functions of the Cartan variables X i by the Schur polynomials where the irreducible representations of the unitary group are labeled by Young diagrams, or equivalently, integer partitions ρ. We provide more details on this in Appendix A.
Observe that the dependence of the Wilson loop on the label α tells us on which of the two supersymmetric cycles the holonomy is evaluated.
By introducing the auxiliary set of power sum variables p s defined as where we introduced two infinite sets of formal time variables τ s,α conjugate to the p s , using the shorthand notation τ α = (τ 1,α , τ 2,α , . . . ). By taking derivatives in times of the generating function and subsequently setting all the times to zero we automatically get all expectation values of the power sum variables and consequently of the Schur polynomials, in other words the WL (α) ρ . Our goal in the following sections will be that of computing recursively such expectation values by making use of matrix models techniques.
For later convenience we also define the exponentiated variables 3 which, as remarked in [10], provide a natural way to describe the generating function as a vector in a representation of two commuting copies of the q-Virasoro algebra (see Section 3.3 for more details on this). As a convenient notational shorthand we will also be using the variable p α with p α = q α t −1 α (not to be confused with the power sum variables {p s }). We remark here that for generic values of the squashing parameters ω α in the region where Im ω 1 ω 2 ≠ 0 (see (B.9)), the contribution of the fundamental antichiral multiplets can be written as .

(2.12)
This contribution can equivalently be obtained (up to a numerical factor) via a redefinition of the CS level and FI parameter together with a shift of the time variables [10]. The corresponding shifts are .

(2.13)
However this transformation becomes singular in the round sphere limit as we discuss in Section 5.
In what follows we will for brevity also use is the integrand of (2.10), where it should be noted that ⟨f ⟩ τ still has a dependence on the times τ 1 and τ 2 (hence the label). With this notation we then have ⟨1⟩ τ ≡ Z (τ 1 , τ 2 ) and more generally In the next section we present the procedure to obtain the q-Virasoro constraints on the generating function Z (τ 1 , τ 2 ) using matrix model techniques.

Derivation of q-Virasoro constraints
For the gauge theory described above, we would now like to derive the q-Virasoro constraints using the trick of inserting a suitably chosen q-difference operator under the integral. This procedure will be very similar to that in [14], the main difference being that now we are working in the logarithmic variables X ∼ ln λ.

Definitions
The finite difference operatorM i,α we will use in order to derive the q-Virasoro constraints is defined as:M where the notation is inspired by [30]. In other words,M i,α corresponds to the operator that sends λ i,α to λ i,α q α (as can be seen from (2.11)). What we now wish to compute is the insertion of under the integral in (2.10) with . . . denoting the integrand and For this reason we will be treating only one copy (i.e. one value of α = 1, 2) at the time. The idea is then on the one hand to compute the action of the operatorM i,α on the integrand, and on the other hand to trade this finite difference operator for a redefinition of the variables X i together with a shift of integration domain (where the details are outlined in Appendix C). We then wish to equate these two computations and obtain the desired constraint.
The motivation for the form of this insertion can be seen as follows. The fraction appearing in (3.2) can be written as where the z-dependence is formal. This enables us to expand the obtained constraint in powers of z, generating a separate equation for each power (i.e. a set of Ward identities in the spirit of those of [31]). Notice here that the summation in (3.4) starts from −1 which in the language of [14] corresponds to considering generic and special constraints simultaneously. This is in analogy with the derivation of the usual Virasoro constraints where one considers the insertion of the differential operator ∂ ∂X i X n+1 i . . . inside of the integral, only now we have to substitute the usual derivative with an appropriate q-difference operator, namely the combination G i,α (λ)M i,α . Furthermore, the precise form of the function G i,α (λ) is necessary in order to introduce the desired pole structure (as shown in Appendix C). We remark that its form is highly reminiscent of that of the Macdonald-Ruijsenaars operator D q,t [11,32] (although the exact relation is yet to be determined).

Computing the insertion
As explained above, what we now wish to evaluate is the following equation where the contour of integration C is taken to be a middle dimensional subspace of C N such that all the integration variables are purely imaginary, i.e. C = (iR) N ⊂ C N (see Appendix C for more details on this).
Here we have introduced the notation (LHS) for the left hand side of the equation and (RHS) for the right hand side so that we can discuss them separately. The (LHS) has been obtained by trading the finite difference operator with a redefinition of the integration variables X i and a shift of the integration domain, which as shown in Appendix C, leaves the integral unchanged. To then evaluate the (RHS) we will compute the action ofM i,α on all the terms above. We remark that while (3.5) holds for physical (real) values of the squashing parameters, all subsequent manipulations of this section are valid for any complex values of the parameters ω 1 and ω 2 .
Starting with the (LHS), (3.5) will hold at each order in z separately, and so we can use the algebraic identity to evaluate the (LHS). Thus using the ⟨. . .⟩ τ notation defined in (2.14). Notice that the first term in the second line is an expectation value of a negative power of λ i,α which we require to be canceled by a similar term on the RHS, as we do not have an interpretation for such terms as differential operators acting on the generating function.
Let us now proceed to the (RHS) of (3.5) by computing the variations, in other words the action ofM i,α , on each of the terms in the insertion and the integrand of (2.10) separately. Starting with the insertion, the variation of We can then use the quasi-periodicity property in (B.7) to compute the variation of the double sineM using which we can evaluate the variation of the measure ∆ Ŝ The variation of the antichirals then becomeŝ where for convenience we introduced P (λ i,α ) as the quadratic polynomial defined by Observe that this polynomial is of degree 2 precisely because we consider the inclusion of two antichiral fields. This fact will be important in Section 4 where we will find that the recursion relates correlators of degree d with those of degree d − 1 and d − 2.
The variation of the CS and FI terms in (2.15) iŝ and finally that of the exponential of the times iŝ We now evaluate the (RHS) in (3.5) by inserting the variations computed in (3.8)-(3.15) above, giving where (following [14]) we have rewritten the expression inside of the average as a sum of residues at the points w = λ −1 i,α , for w an auxiliary complex variable. Now we use the fact that w is a point on the Riemann sphere to move the contour in such a way that it encircles the poles at w = 0 and w = ∞, 4 (with opposite orientation) instead of the poles at w = λ −1 i,α .
The (RHS) can thus be rewritten as where we are now able to bring the matrix model average inside of the w integral.
First we compute the residue at w = ∞ in (3.17) by substituting F (w) with its power series expansion so that we obtain In order to determine the coefficients F n we need to specify the value of the CS level κ 2 .
For the result to be non-vanishing we first of all require κ 2 ≤ 2. Secondly, the choice of κ 2 has to be such that the residue at w = ∞ yields a term which can cancel the expectation value of λ −1 i,α appearing in (3.7). The only consistent such choice (that does not introduce any higher negative powers of λ i,α ) is which will be used in what follows. Being in a neighborhood of w = ∞ we can use the identity to rewrite the last term in F (w) and we immediately see that all coefficients F n vanish for n > 0 so that the only contributions to (3.19) are those for n = −1, 0. An explicit computation gives We now consider the other residue in (3.17), namely the residue at w = 0. We first rewrite the integral as where F n are now the coefficients of the power series expansion of F (w) around w = 0 and remainder is the part of the series which only contains negative powers of z of degree less than −1. Given that we are only interested in the q-Virasoro constraints for n ≥ −1, such spurious terms can be neglected in the derivation of the main equation and therefore we will not be interested in writing their particular expression.
Next, as we are working in a small enough neighborhood of w = 0, we can use the identity to rewrite the residue as Finally, we can combine both residues and plug everything back into the original equation (3.5), to obtain (3. 26) In order to have the cancellation of the term ⟨1 λ i,α ⟩ τ as discussed above, we set the value of κ 1 accordingly 5 27) which leads to the constraint equation Thus, what we did is to rewrite the finite difference equation (3.5) as a differential equation in the time variables τ s,α for the generating function Z(τ 1 , τ 2 ). Upon expanding this equation in powers of z (for n ≥ −1) we obtain a set of differential constraints which we can interpret as an explicit representation for the q-Virasoro algebra (see Section 3.3). From now on we will refer to (3.28) as the (combined) q-Virasoro constraints.
In Section 4 we will provide a recursive solution for this set of constraints.

Free field representation
Before attempting to solve equation (3.28) we want to provide an algebraic description of the constraints and their relation to the representation theory of the q-Virasoro algebra as previously studied in [10,11]. The generating function defined in Section 2 can be interpreted as a highest weight vector in a module over two commuting copies of the q-Virasoro algebra, one for each value of the label α. Each copy of the algebra is generated by the operatorsT n,α for n ∈ Z with which we can express the q-Virasoro constraints aŝ T n,α Z(τ 1 , τ 2 ) = 0, n ≥ 1 (3.29) expressing the condition that Z(τ 1 , τ 2 ) is indeed a highest weight vector annihilated by all the positive generators (the generatorsT 0,α are diagonal on the highest weight and they give simple eigenvalue equations when acting on the generating function).
Notice that here we use the "hatted" notationT α (z) instead of the standard non-hatted current of [11] to stress that the representation of the algebra is deformed by the introduction of the antichiral fundamental multiplets (2.12), which amounts to the time shifts (2.13).
It is then customary to package the full set of generators into a stress tensor currentT α (z) asT α (z) ∶= n∈ZT n,α z n (3.30) 5 This can be compared to the value of κ1 in [10] where the additional terms − ω 2 + µ+ν 2 are generated by the inclusion of the fundamental chiral multiplets and we set their parameter α = 0 (not to be confused with our index α). so that the constraints can be collectively rewritten aŝ where Pol α (z) is a function whose power series expansion only contains non-positive powers of z. By expanding in powers of z on both sides of the equation, one recovers the action of each of the generators of the algebra.
What we want to do now is to interpret equation (3.28) as a concrete representation for the algebraic identity (3.31). In order to do that, we first introduce the following representation for the Heisenberg oscillators of [11, Section 4] satisfying the algebra [a n,α , a m,α ′ ] = n By using this explicit free field representation and also introducing the function ψ α (z) as we are able to rewrite (3.28) as 35) where the currentT α (z) takes the form of the differential operator and remainder is identified with the part of ψ α (z)Pol α (z) with powers of z of degree less than −1.
Comparing (3.36) with the formula for the current of [11] we observe that the only difference is the multiplicative factor P (q α z) appearing in front of the second term. This deformation is due to the presence of the two anti-fundamental flavors of masses µ and ν on which the polynomial P depends through the coefficients A, B. As a direct consequence we observe that the constraint equation forT −1,α is also modified aŝ (3.37) The equation forT 0,α instead does not depend on the deformation and gives the usual eigenvalue equationT 0,α Z (τ 1 , τ 2 ) = p that we expect from a highest weight module representation.

Recursive solution
The goal of this section is to solve the q-Virasoro constraints in (3.28), where by solve we mean to recursively determine all the normalized correlators C 1 ... n;k1...km of this model. The correlators are defined as We assume that the partition function admits the formal power series expansion: where Z d 1 ,d 2 (τ 1 , τ 2 ) has degree d 1 and d 2 with respect to the operators ∑ ∞ ∂ ∂τ d 2 ,2 , respectively. We then use the definition (A.9) of the symmetric Schur polynomial s {m} (p 1 , . . . , p m ) for a given symmetric partition {m} in order to extract the coefficient of z m , m = −1, 0, 1, . . . in (3.28). When doing so, we only consider terms of a particular degree d 1 in times τ ,1 and degree d 2 in times τ k,2 .
By inserting formula (4.2) into the q-Virasoro constraint (3.28) and choosing α = 1 for definiteness, we get The corresponding equation for α = 2 is completely analogous but it is d 2 that is shifted by m.
Given any two partitions ρ = {ρ 1 , . . . , ρ • } and σ = {σ 1 , . . . , σ ⋆ }, where ρ • and σ ⋆ indicate the last components, we wish to compute the correlator C ρ;σ ≡ C ρ 1 ...ρ•;σ 1 ...σ⋆ , with the identifications m + 2 = ρ • , d 1 + 2 = ρ and d 2 = σ . To extract the correlator we are interested in, we apply the operator to (4.3), namely we differentiate with respect to the corresponding combination of times and then set all of them to zero. Finally, we obtain which can be used to determine the correlator C ρ;σ in terms of the lower order correlators on the right hand side. Here γ and η are partitions and we used the formulas of Appendix A to rewrite explicitly the symmetric Schur polynomials s {m} ({p k }) in (4.2). The sum over subsets η ⊆ ρ ∖ ρ • means that we are summing over all the sub-partitions η of ρ not containing its last part ρ • . An analogous equation holds for simplifying the multi-index σ of the correlator.
Every correlator is uniquely determined by repeated application of (4.5). Therefore, if a non-trivial solution to the over-determined system of q-Virasoro constraints (3.28) exists, then it must be given by (4.5). While the proof of consistency of (3.28) is out of the scope of this paper, we did check (up to degree 6) that it is indeed consistent.
The recursion in (4.5) treats the two copies α = 1, 2 separately. From this observation, one might be lead to believe that the correlators factorize, although this is not obvious from the form of the partition function (2.7). We now prove that this is indeed the case. First of all suppose that all correlators for partitions up to a certain order d factorize into the product of correlators for α = 1 and α = 2 as Then, if we consider a partition ρ of order d + 1, all correlators in the right hand side of the recursion formula (4.5) factorize as per our assumption and therefore we can collect a common factor of C ∅;σ . This immediately implies that the ratio C ρ;σ C ∅;σ can be written entirely in terms of correlators of the form C ρ ′ ;∅ for ρ ′ some partition of order at most d. In other words, the correlator on the left hand side (which is of order d + 1) must also factorize as prescribed in (4.6). Finally we need to show that the initial data of the recursion factorizes as well. However, the only correlator needed to fix this data is the trivial correlator C ∅;∅ whose value can be set to 1 as a choice of overall normalization which corresponds to looking at the normalized correlators. The factorization of all correlators then follows by induction, as well as that of all Wilson loop expectation values for any pair of partitions ρ, σ. While non-obvious from first principles, this factorization property is similar to the holomorphic blocks factorization of [33,34].
With the help of the recursion we find the first few correlators , (4.9) where we recall that A and B are the coefficients of the polynomial P (λ) given in (3.13).
For the first few correlators of Schur polynomials we manifestly get , , which can be translated to expectation values of Wilson loop operators via the relation Computation of any other Wilson loop expectation value goes through in exactly the same way.
In the next section we study several limits of the recursive solution outlined above.

Limits
There are several interesting limits that can be taken starting from the result obtained in the previous section.

Round sphere
Firstly, one can consider the round sphere limit which corresponds to letting ω 1 = ω 2 , or equivalently q α → 1 for both values of α. In this case the value of t = t 1 = t 2 = e 2πiMa is kept fixed and the q-difference equation (3.28) takes the simpler form: 1) where, since there is no distinction between the two copies α = 1 and α = 2, we can define a single set of time variablesτ s = τ s,1 + τ s,2 . The partition function Z S 3 (τ ) of the gauge theory on the round sphere is a solution of this equation. Observe also that in this case the two copies of the q-Virasoro algebra collapse to the same giving rise to just one set of constraints.
Furthermore we note that in this limit the shift in the time variables, which can be used to reproduce the contribution of the fundamental antichiral multiplets, becomes singular (see (2.13)). This is also evident from the fact that in the round sphere limit we have Im(ω 1 ω 2 ) = 0 which is the region where the double sine cannot be factorized into a product of q-shifted factorials (see (B.9)) and therefore (2.12) does not hold.
Besides, in this limit the current in (3.36) is such that it only contains derivatives in times. This can be seen from the observation that only half of the Heisenberg oscillators in (3.32) will remain, namely the positive modes. This renders the algebra of both the Heisenberg oscillators and the q-Virasoro generators abelian. We remark that this also corresponds to the Frenkel-Reshetikhin limit W 1,t of [35].

Gaussian matrix model
Secondly, by making the specific choice of masses such that v α = −u α = (1 − q α ) −1 2 , the contributions of the fundamental antichirals in (2.12) simplify to q-Pochhammer symbols using (B.9). For instance, for α = 1 where the ambiguity in the sign originates from interpreting the minus sign as an exponent. It can be noted that this depends only on ω 1 and not ω 2 , and so we expect that this particular choice can only be made for one value of α at a time. By setting the masses as in (5.2) the contribution of the fundamentals introduces a q-exponent factor in the first set of variables λ i,1 of the form which reproduces the potential of the (q, t)-Gaussian 6 model studied in [14]. The dependence on the second set of variables λ i,2 however, cannot be put in the form of a q-exponent function for the same values of masses.

U (1) gauge theory
Thirdly, one can consider the N = 1 limit of the localized partition function on S 3 b , which is the case originally studied in [33]. The theory becomes that of a single U (1) gauge group with two (anti)fundamental flavors while the adjoint multiplet becomes a singlet and decouples. In this case the integral in the partition function becomes 1-dimensional and the measure ∆ S (X) in (2.7) is simply 1. As a consequence of this simplification, the correlators in this limit do not depend on the adjoint chiral mass M a and therefore t α does not appear in the formulas. As is expected for a rank 1 matrix model, all correlators associated to partitions of the same degree become equal to each other. For example we have

(5.4)
Moreover, one finds that all Schur polynomials vanish identically except for those associated to completely symmetric representations, i.e. to those partitions of length 1, so that we have for all m ≥ 0 where m is the U (1) charge of the representation.

Conclusion
In this paper we developed an explicit algorithm that allows one to calculate supersymmetric Wilson loop averages in N = 2 YM-CS theory coupled to specific matter on a squashed sphere S 3 b for any value of b, including b = 1. The heart of the algorithm is the Ward identities (3.28) for the corresponding matrix model, which produce the recursive equation (4.5) that expresses power sum monomial correlators through simpler ones.
Due to technical reasons we restrict our attention to specific values of CS level κ 2 and FI parameter κ 1 . Our current understanding suggests that such technical limitations reflect an underlying deeper conceptual difficulty related to the possibility to solve the matrix model only around a specific "expansion point", similarly to how the usual Hermitian matrix model is completely solvable by expanding around a Gaussian potential. Extending the present result to generic values of CS level and FI parameter is certainly an interesting question which requires further investigation.
The present analysis is the natural continuation of the works [10] and [14]. In principle the present work could be extended to 3d N = 2 unitary quiver gauge theories (which include the ABJ(M) model [36,37] and its deformations) as described formally in [10] with the Ward identities related to W q,t (Γ) algebras of [38]. The main problem is how to disentangle the algebraic and analytical issues in this generalization. Within the formal free field theory representation all quiver theories can be treated uniformly. However we do not know how to analyse q-difference operators, poles and contours in a uniform way. At the moment we can only do it model by model and it is a quite time-consuming analysis.
In a larger context it would be nice to understand how to solve q-Virasoro constraints in a more general situation. The generating functions for 3d and 5d gauge theories are naturally related to q-Virasoro constraints (and its relevant deformations [39]). For example, the Nekrasov generating function can be thought of as "N = ∞" (q, t)-deformed matrix models, see the Kimura-Pestun representation of Nekrasov function [38].
Last but not least, while from the point of view of U (N ) gauge theory Schur polynomials are the objects to consider, from a purely matrix model point of view a different kind of polynomials -the Macdonald polynomials -are much more natural. They are the orthogonal polynomials for the (q, t)-Vandermonde measure, and their average is actually a very simple factorized expression [40,41] -a property known as "character expansion" [41]. Therefore, understanding the gauge theory significance of Macdonald polynomials is an important problem for future research.

A Schur polynomials and partitions
Let us begin by introducing the notation for integer partitions. Consider the partition γ = {γ 1 , . . . , γ • } where γ 1 ≥ ⋅ ⋅ ⋅ ≥ γ • > 0 are positive integer numbers. Then γ is the degree of the partition γ γ ∶= where on the right hand side we sum over all partitions γ and the variables p k are usually identified with traces of powers of some N × N matrices Φ In the main part of the paper we use the identification where α = 1, 2 labels the two sets of variables (2.11) of the matrix model.
In the particular case when one adopts the plethystic substitution τ k = z k , then on the right hand side of (A.5) the summation collapses to the symmetric partitions γ = {m} only, i.e. those of length 1. Using the formula for symmetric Schur polynomials More generally, Schur polynomials for arbitrary partitions γ can be expressed using determinants as The first few Schur polynomials are given by (A.11) We end this section by mentioning that integer partitions and Young diagrams have a deep relation with the representation theory of the symmetric group S n and, more importantly, the representation theory of U (N ) (and SU (N )). There is a 1-to-1 correspondence between partitions and irreducible representations 7 of the unitary group and because of this it follows that Schur polynomials are exactly the irreducible characters of U (N ). For a given irreducible representation R γ associated to the partition γ and a group element Φ, we have

B Special functions
Here we collect a few well known facts about the special functions that we employ throughout the paper. For additional details on the multiple sine functions we refer the reader to [43].
We first introduce the q-Pochhammer symbol, or q-shifted factorial, which is defined as with x ∈ C and q < 1. The analytic continuation to the region q > 1 is given by the formula Secondly we introduce the double sine function. For ω ≡ (ω 1 , ω 2 ) ∈ C 2 with Re(ω α ) > 0 and z ∈ C, the double sine function is defined by the regularized product where ω = ω 1 + ω 2 . Let Λ be the semi-lattice then the double sine is a meromorphic function with zeroes and poles at zeroes ∶ z = −Λ , as shown in Figure 1.
It can be shown that it satisfies the reflection identity as well as the first order finite difference equations The highest weights γ of the irreducibles of U (N ) are those that are integral (γi ∈ Z) and dominant (γ1 ≥ γ2 ≥ ⋅ ⋅ ⋅ ≥ γ N ). Moreover the double sine function is related to the hyperbolic gamma [24] function Γ h by S 2 (z ω) = Γ h (z; ω 1 , ω 2 ) −1 . (B.8) Additionally, for Im(ω 1 ω 2 ) ≠ 0 one can express the double sine function as a product of two q-shifted factorials [44] S 2 (z ω) = e

C Difference operator and shift of countour
Let us consider the integral in (2.10). If we assume that the squashing parameters ω 1 , ω 2 are both real and positive, then we can fix the integration contour C to be a middle dimensional imaginary contour inside of C N . More precisely let (iR) N be the subspace defined as (iR) N ∶= X ∈ C N Re(X i ) = 0, for i = 1, . . . , N ⊂ C N , (C.1) then the integration contour is prescribed by supersymmetric localization to be C = (iR) N . With this choice the integrand falls off fast enough at complex infinity in the right-halfplane and one can compute the integral using the residue theorem by closing C to the right. If there are poles lying exactly on the imaginary axis (of one of the variables) we give them a small real part or equivalently we deform slightly the contour to their left.
For n ≥ −1, we want to show the following identity between integrals: where we consider one index i and one order of z at a time in (3.5) (using (3.4)). The nontrivial part of the identity is the second step where the shifted contourM −1 i,α C is deformed to the original contour C. The result of the integration will then be unchanged only if there are no poles in between the two contours. We now show that this is in fact true.
In order to show that there are no poles between C andM −1 i,α C, we first notice that the shift operatorM i,α only acts on the variable λ i,α while leaving all other λ j,α for j ≠ i unchanged. This means that we can simplify the problem by just focusing on the integration in the complex plane of X i . Having fixed the index i we can regard the other variables as background imaginary parameters. For definiteness we also take α = 1.
The poles that we should worry about are those coming from the function G i,1 (λ), the measure ∆ S (X) and those coming from the antichiral multiplets. The latter are the poles of S 2 (−X i + µ ω) −1 and are situated at the positions which we assume to be to the right of the integration contour. 8 After the shift, the poles of the functionM i,1 S 2 (−X i + µ ω) −1 are situated at so that clearly they are all still to the right of the contour C. A similar analysis follows for the mass ν.
Next we look at the Vandermonde measure ∆ S (X). Observe that the term ∏ j≠i S 2 (X i − X j ω)S 2 (X j − X i ω) in the numerator only has zeros and no poles because of the identity S 2 (z ω)S 2 (−z ω) = −4 sin πz ω 1 sin πz ω 2 . (C.5) The zeroes at X i = X j + ω 1 Z then cancel exactly the poles of the function G i,1 (λ) coming from the denominator in (3.3).
The term ∏ j≠i S 2 (X i − X j + M a ω)S 2 (X j − X i + M a ω) in the denominator of ∆ S (X) is a bit more involved. When one of the X j is infinite then the part of the Vandermonde which depends on X i and X j goes to 1 and there are no poles because the adjoint mass becomes negligible and the numerator cancels with the denominator (in this regime the function G i,1 (λ) is subject to the same kind of cancellation between poles and zeroes). When the X j are finite we can regard them as finite imaginary shifts in the poles where for simplicity we assume that the adjoint mass M a has a small and positive real part so that the contour separates the poles at X j + M a + Λ (on the right) from those at X j − M a − Λ (on the left). The insertion of the function G i,1 (λ) has the effect of removing the poles situated at the boundary of the semi-lattice X j − M a − Λ. If we now apply the shift operatorM i,1 we get that the shifted positions of the poles are: We observe that those poles that would have crossed the imaginary axis because of the shift are precisely those at X i = X j − M a − ω 2 Z ≥0 , i.e. those that are canceled by the numerator of G i,1 (λ), so that the contributions to the residue computation are unchanged by application ofM i,1 . In conclusion we have that the contourM −1 i,1 C can be deformed to C without crossing any poles and therefore the identity (C.2) holds.