Solving q-Virasoro constraints

We show how q-Virasoro constraints can be derived for a large class of (q, t)-deformed eigenvalue matrix models by an elementary trick of inserting certain q-difference operators under the integral, in complete analogy with full-derivative insertions for β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document}-ensembles. From free field point of view, the models considered have zero momentum of the highest weight, which leads to an extra constraint T-1Z=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{-1} \mathcal {Z} = 0$$\end{document}. We then show how to solve these q-Virasoro constraints recursively and comment on the possible applications for gauge theories, for instance calculation of (supersymmetric) Wilson loop averages in gauge theories on D2×S1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D^2 \times S^1$$\end{document} and S3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S^3$$\end{document}.


Introduction
Over the last 40 years, matrix models have attracted enormous attention in mathematical physics. Besides their direct applications, matrix models provide a useful playground for quantum field theory and as such they have been extensively studied (see [1] and references therein).
The simplest examples of matrix models are ordinary integrals over the space of finite dimensional matrices. However, the real interest is in more complicated, deformed, examples, where the link to the integrals over matrices is less direct or even completely lost. The central question is then: which (hidden) structures do survive the deformation [2]? The recently renewed interest in deformed matrix models is fueled by applications in gauge theories (calculations using localization technique [3] often lead to an effective description in terms of matrix models), as well as the theory of symmetric polynomials and quantum groups [4,5]. The deformation that is immediately relevant to supersymmetric gauge theories is the so-called (q, t)-deformation. In this paper we show that certain standard tools available in the non-deformed cases can also, in fact, be used in the (q, t)-deformed case.
The central object in any matrix model is the partition function Z, which is, usually but not always, a formal power series in infinitely many "time" variables {t k }. Coefficients of Z with respect to {t k } are called correlators and one of the questions for each given class of models is how to find these correlators, ideally in a fast and efficient way.
In the case of non-deformed models, exemplified by the Hermitian 1-matrix model [1,Section 3.2], there are many answers to the question of how to find correlators. One way is to use the character decomposition formulas, which are also conjecturally available in the deformed case [6]. Another method is to derive Ward identities and use them as equations for correlators. Yet another way is to derive the so-called loop equations from the Ward identities and further recast them into the algebrogeometric form known as the spectral curve topological recursion [7][8][9][10]. All these different ways, which can be thought of as manifestations of different hidden structures of the model, are straightforwardly linked to one another in the non-deformed case (and in the case of Hermitian Gaussian 1-matrix model it is especially easy to trace all the connections). However, after (q, t)-deformation, each of the hidden structures changes in its own peculiar way. As a result, relations become less apparent and so each structure needs to be studied by itself.
Currently, some of the deformed structures are understood more and some are understood less. Character expansion is confirmed by computer experiments, yet its derivation from first principles is a challenging open problem [6]. The (q, t)generalization of the topological recursion is completely unknown (but, on the other hand, topological recursion for β-ensembles is well understood, see [11]). It is not even clear how to distinguish contributions of different genera in a correlator. The situation with the (q, t)-deformed Ward identities, which are the main focus of this paper, is the following. In the non-deformed case, there are two different ways to derive the Ward identities. The first one, which is more conceptual as it exhibits connection between matrix models and conformal field theories, relies on representing the partition function of a given model as an average in some free field theory [12]. The second one, which looks more like a clever trick, consists of insertion of a full derivative under the integral [1, Section 2.1]. The (q, t)-version of the first derivation is well known [12]. From it, we know that for a large class of (q, t)-deformed models Ward identities form the so-called q-Virasoro algebra [13]. However, the derivation of the q-Virasoro constraints using the insertion of a full derivative (which should be substituted with a full q-difference operator) is not known in general, but there are some exceptions. In [14], it is done in the case of (q, t)-deformed Selberg integral. The (q, t)-deformed Selberg integral is a very special point in the space of all (q, t)-deformed matrix models as it is the only matrix model for which the character expansion is proved [4, Chapter VI, Section 9, Example 3]. The sketch of the generic derivation can be found in [15,Section 2.1], though the details are not spelled out. In this paper, we do derive the q-Virasoro constraints using the insertion of a full difference operator (details are in Sect. 4) for a large class of models (defined in Sect. 2).
As it often happens, working out the details brings interesting surprises. In order for the equations to be sufficient to determine all correlators in the simplest case of Gaussian (q, t)-deformed model, we need to consider one "additional" equation (see Sect. 4.2.2), which is the (q, t)-analogue of the string equation [9,Introduction]. When we spell out what this equation means in the language of free fields, it turns out to be nothing but the condition T −1 Z = 0 (usually the partition function is annihilated only by positive q-Virasoro generators, see Sect. 3). This is an extra equation that is valid only for partition functions with zero vacuum momentum.
Having derived q-Virasoro constraints we proceed to solving them. Namely, we describe, how to use them to find all correlators of the model. In the non-deformed case, thanks to the simple form of the Virasoro constraints (they are second-order differential operators) it is obvious how to convert them into efficient recursive procedure, that allows to calculate each given correlator in polynomial time. The form of q-Virasoro constraints is more complicated, and it is far less obvious how to use them to get an efficient recursion. Still, we find such recursion in Sect. 5 (see Eq. 5.4).
This efficient procedure for evaluation of correlators is of immediate interest for gauge theory applications. In gauge theories, one is interested in supersymmetric Wilson loop averages. On the (q, t)-matrix model side, these averages correspond to the correlators of Schur functions (even though the Macdonald functions are more natural from the (q, t) point of view). Each Schur function is a finite linear combination of power sum monomials; therefore, one can use the recursion procedure of Sect. 5, to evaluate any supersymmetric Wilson loop in any gauge theory, that corresponds to some matrix model from the class we define in Sect. 2.
Last but not least, from the free field derivation, one can see that (q, t)-deformed models have a vast hidden freedom which is not present in the non-deformed case: one can insert arbitrary q-constants under the integral. In our derivation using difference operators we also see that Ward identities for discrete and continuous Gaussian matrix models [Eqs. (2.6) and (2.7), at the special choice of masses and Fayet-Iliopoulos parameter, respectively] are the same (and hence the correlators are the same). From the gauge theory point of view, this looks like a non-trivial duality between the corresponding gauge theories. So, the q-Virasoro viewpoint is potentially very fruitful in the search for dualities.
The rest of the paper is organized as follows. In Sect. 2, we define the class of models we study and provide the two models we use as examples. In Sect. 3, we give the matrix model background and the free field realization of the q-Virasoro algebra, and in Sect. 4 we move on to the derivation of the constraint equations which the models have to satisfy. In Sect. 5, we consider the recursive solution to correlators and in Sect. 6 we see how this can be applied to gauge theory. Finally in Sect. 7, we conclude and suggest further interesting directions of investigation. Definitions of special functions and details of computations are left to the appendices.

Definition of the class of models
Let us consider a class of models of dimension N where the partition function is given by for some appropriately chosen contours {C i }, with integrand Here (x; q) ∞ = ∞ k=0 (1 − xq k ) is the q-Pochhammer symbol and c q is of the form where α ∈ C [later to be interpreted as momentum of CFT states, see (3.12)] and λ q (x) is a q-constant. 1 The parameters t, q, β ∈ C are related via t = q β , while This is a natural (q, t)-generalization of eigenvalue matrix models and of βensembles: indeed, the first factor in (2.2) plays the role of (q, t)-deformed Vandermonde measure, the second introduces time variables {t k } and the third encodes the common potential for the "eigenvalues" x i . In gauge theory applications, N f is the number of (anti-)fundamental masses and the parameters of the potential m f are the respective fugacities. Note that the product Of course, this is not the first time this class of models makes an appearance in the literature. For example, in pure mathematics they appeared in the context of constant term evaluations and Macdonald conjectures (see [16] for a review). In the physical context they played an important role in the description of AGT dualities [17], trialities between conformal field theories and gauge theories in three and five dimensions [18,19] and more recently hidden symmetries of network-type matrix models and their relation to DIM algebras [15].
Within this class of models there are two examples that we will frequently use for illustration purposes. The first one is the discrete (q, t)-Gaussian model, obtained by taking N f = 2 with two opposite masses ±m having the special values ±m = ±q(1 − q) 1 2 . In the limit q → 1 the product of the two q-Pochhammer symbols becomes e − x 2 2 and the standard Gaussian potential is recovered. This justifies the name (q, t)-Gaussian model as the simplest deformation of the Gaussian 2 model. This model was axiomatically defined in [6] via its Macdonald averages, but here we give the explicit integral representation. The function c q (x i ) for the (q, t)-Gaussian model can be written as is a particular q-constant. Alternatively, the same model can be expressed in a more concise form using the notion of a Jackson q-integral-a discrete q-analogue (A.5) of an integral-via summing over residues (A.9), provided that contours are chosen as in Fig. 1: These two descriptions-continuous and discrete-provide two useful and complementary perspectives on one and the same model. The second example that we consider is the 3d N = 2 gauge theory partition function on D 2 × S 1 where the non-trivial decomposition was found in [20] and the q-Virasoro interpretation was found in [18] and extended in [12]. The holomorphic block integral (omitting factors of 2πi) for this gauge theory takes the form where the integration is over middle dimensional cycles c in (C × ) rkG , m a is the adjoint fugacity, m f and m f are the fugacities for the fundamental and anti-fundamental and κ 1 is the Fayet-Iliopoulos parameter. As given in [12], we have the identification κ 1 = √ β(α 0 + √ β N − Q β ) (for some initial momenta α 0 ∈ C) between the gauge theory and the q-Virasoro side, and the model is then reproduced by with the generating function We will now move on to review the properties of the free field realization of q-Virasoro algebra that will be needed to understand the above constraints.

(q, t)-deformed Virasoro matrix model and free fields
In this section, we recall the necessary facts about the (q, t)-deformed eigenvalue models and their free field realization, where details can be found in [1,12,21,22]. It can be noted that any model of the form (2.1) can be generated using the free field representation given below, and so there is no discrepancy between the different point of views.
The q-Virasoro algebra is given by the set of generators {T n , n ∈ Z}, subject to relations or equivalently, in terms of generating currents where and the parameters p, q, t ∈ C and p = q/t. The free field representation of the q-Virasoro algebra is given in terms of the following Heisenberg oscillators using which the generator current takes form (3.5) Explicit expressions for the generators are given by Here we let and is the Schur polynomial in symmetric representation [n], defined by (A.4) and in particular s 0 = 1 and s 1 (A 1 ) = A 1 . We will also use the time representation of the Heisenberg algebra in terms of time variables {t k }: We can then define the screening current x −n (q n/2 −q −n/2 ) a n : e √ βQ x √ βP , (3.10) which has the nice property of commuting with the generators up to a total q-derivative for some operator O n (x). Using this screening current we can then consider the state for some momentum α ∈ C where the integration contours are chosen such that upon the insertion of a total derivative D, the integral vanishes: Here, it can be mentioned that this requirement can be shown to hold for a wider class of models which can have measures other than N i=1 dx i x i , e.g., one can use either the Jackson integral or integral over Since the vacuum |α is annihilated by positive q-Virasoro generators, the state (partition function) Z{t k } in Eq. (3.12) also satisfies the q-Virasoro constraints for some function f (z) holomorphic at z → 0. Since the vacuum is an eigenvector of T 0 [13] the partition function is also an eigenvector of T 0 with the same eigenvalue, Additionally, restricting to the subset of states which has zero momentum, α = 0, we have an additional equation (verified in Appendix E) Two points deserve to be mentioned here. First, from the point of view of full qderivative insertion, the additional equation T −1 Z{t k } = 0 arises quite naturally (see Sect. 4.2.2). The requirement of zero momentum there takes the form of the vanishing of the coefficient in front of some complicated term, which is otherwise unmanageable. Once this term is gone, at least for (q, t)-Gaussian matrix model, we can recursively determine all the correlators.
Second, from the point of view of the q-Virasoro representation theory, the condition T −1 Z = 0 (which can be rewritten as [T 1 , T −1 ]Z = 0) is necessary to have the so-called singular vector at level 1. In a sense, it is the simplest irreducible representation of q-Virasoro. It would be interesting to know, which matrix models (in a sense of explicit integral representation) correspond to more complicated q-Virasoro representations.

Ward identities via q-difference operator insertion
We now proceed to the derivation of the Ward identities. First, in Sect. 4.1, we define the q-difference operator that is to be inserted under the integral. Then, in Sect. 4.2, we gradually rewrite the q-difference insertion as the action of shift operator in times {t k }, which culminates in Eq. (4.20)-the main result of this section. Finally, we compare the resulting Eq. (4.20) with the q-Virasoro equations (3.14) confirming that they match. On the way we discover that in order to proceed we need to impose certain conditions on the parameters of the model, in particular on the function c q (x) and on the integration contours for the eigenvalues.

The q-difference operator
We want Ward identities to be a corollary of the fact that insertion of the suitably chosen full q-difference operator under the integral vanishes, i.e., that where the q-difference operator is and F(x) as given in (2.2). Here, the dependence of G i (x) on z is understood formallywe will obtain an equation for each coefficient in front of particular positive power of z. For Eq. (4.1) to hold, we must impose certain restrictions on the integration contours and the potential term c q (x i ). Namely, let us consider terms in (4.1) that contain operator x † . This operator can be traded for the change of the integration contour as follows Assuming that the complex parameter q satisfies |q| < 1, the contour C i /q will encircle C i as shown in Fig. 2. So, in order for (4.1) to hold, we must be able to shrink the contour C i /q back to C i so that we can take the sum over i inside again.
We consider the implications of this condition in the following subsections.

From full derivative to shift operators in times
As we now have ensured that the insertions we consider are well defined and vanishing, we will proceed with rewriting them as some difference operators acting on times {t k }, in complete analogy with the non-(q, t)-deformed case.

The generic constraint
The integrand of the left hand side of (4.5) can be rewritten as (see Appendix B) After we integrate over the eigenvalues, the left hand side of (4.6) can be rewritten as provided we can shrink contour C i /q back to C i . We can do this if there are no poles of the subintegral expression in the shaded region. Let us look at different parts of the subintegral expression in detail.
As we said before, we treat the first part of G i (x), 1 1−zx i , as a formal power series in non-negative powers of z, therefore, it does not contribute any poles. The second can possibly have poles at x j = x i for j = i which could appear in the shaded area. However, the q-Pochhammer in the numerator of F(x), k =l (x k /x l ; q) ∞ will precisely cancel these poles, and so there is no additional contribution from G i (x) as we shrink the contour.
Then, let us consider F(x). The first part, k =l can have poles for noninteger β. However, these poles are canceled by the numerator of G i (x). In the second part of F(x), the function c q (x i ) can, in principle, have poles in the shaded region. One therefore needs to show that this does not happen for the model one is interested in, otherwise Eq. (4.1) would acquire a non-trivial right hand side.
For the (q, t)-Gaussian model, the function c q (x i ) is given by (2.4) and does not pick up any poles as we shrink the contour from encircling (−ν/q, ν/q) to encircling (−ν, ν): indeed, the only possible pole in this area is the point x = ν/q (a pole of f λ ) but it coincides with a zero of the q-Gaussian measure (x 2 q 2 ν −2 ; q 2 ) ∞ and therefore is just a regular point.
For the gauge theory on D 2 ×S 1 model, the argument is analogous: for a given fugacity m f as we shrink the contour from encircling Having demonstrated that for these two models there are no contributions from poles in the shaded area, from now on we proceed assuming this is the case, so that We then use the following algebraic identity (verified in Appendix C) Let us now move to the right hand side of Eq. (4.6) where we also take an integral over eigenvalues. We want to interchange the integral over eigenvalues with operation of taking the residue. But we cannot do this, since the points we take residues at manifestly depend on the eigenvalues. Therefore, we change the sum over residues at ω = q/x i for the sum of residues at ω = 0 and ω = ∞ (details are in Appendix D). Then, we can safely exchange integral over eigenvalues with the residue operator.
We can only do the above transformation provided that the subintegral expression has no other poles except q/x i , 0 and ∞. If there would be extra poles they would contribute to the right hand side of the Ward identities. Then, Ward identities would no longer be q-Virasoro constraints, at least not in the representation (3.9). So, we have an additional requirement that c q (1/ω) c q (q/ω) should have no poles except at 0 and ∞. This requirement should be checked for every concrete model one is interested in. For instance, for the (q, t)-Gaussian matrix model, it is straightforward to see that this is, indeed, the case. For the gauge theory on D 2 × S 1 , however, the fugacity for the fundamental, m f , that appears in the denominator of (2.7) could cause such poles. Therefore, the dependence on the fundamental fugacities can be studied in the framework of q-Virasoro constraints only as a formal power series. Luckily, the dependence on the anti-fundamental fugacitiesm f need not be understood formally, one can study it non-perturbatively.
Hence, Eq. (4.6) becomes (after some massaging that involves calculating residue at ω = ∞) (4.11) Using the explicit form of the generating function given in (2.1) and the simple fact that 3), we can write the above in terms of the partition function Z using shifts in the time variables {t k }: which is then our first constraint equation.
Here, we have rewritten certain products as exponents of logarithms, for instance Let us proceed to the second constraint.

The special additional constraint
There is another insertion of full q-derivative that one can consider that leads to an additional equation for the partition function. In the case of (q, t)-deformed matrix model, this equation is necessary to make the system of equations closed and to be able to find all the correlators. In the classical limit this equation goes into the string equation. To obtain it, we begin with the relation (verified in Appendix B) where nowG Proceeding similarly to above, the first step is to take a contour integral of the left hand side and use the algebraic identity We then transform the right hand side by changing the residues from ω = q/x i to 0 and ∞ and swapping integration over eigenvalues with taking the residue. Writing the result using Z{t k } we get after some calculation Note that the complicated first term vanishes in the case α = 0; this will be useful later. Otherwise, it is not clear what to do with it: it is an average of trX −1 and as such it is not expressible as any differential operator in {t k } acting on Z. Perhaps, the set of times could be enlarged to be able to also generate correlators of traces of negative powers of X ; however, investigation of this possibility is beyond the scope of the present paper.

Combining generic and special constraints
We then want to combine our two constraint equations. The key point here is observation that the first constraint in Eq. (4.12) can be written as which the generating function Z{t k } must satisfy.
As an example and for comparison with earlier work on q-Virasoro, let us consider the special case with no potential, i.e., N f = 0 and α = 0. It can be seen that in this case the above constraint takes a simple form for some coefficients c m . Comparing this to the construction in [13] and the q-Virasoro constraint in Eq. (3.14) (with z → 1/z), we are considering with c 0 and c 1 being the coefficients of the z 0 and z −1 terms, respectively, and To then extract the action of T 0 and T −1 on Z{t k } we expand the generator current T (1/z) = n∈Z T n z n and compute from which we can identify the eigenvalue equations This can also be obtained using the free field representation of [12] as shown in Appendix E, and so these two viewpoints are equivalent.

Recursive solution
The aim of this section is to show how q-Virasoro constraints can be used to recursively determine all normalized correlators of a given eigenvalue model. We will use the (q, t)-Gaussian matrix model as an example, but models with more involved potentials can be treated similarly-they just require more initial conditions to start the recursion.
This situation is completely analogous to the usual, non-(q, t)-deformed case, where one needs to choose a set of integration contours (i.e., the "phase" of the model) to define its correlators [23]. The q-Virasoro constraints (4.20), specialized for (q, t)-Gaussian matrix model, read where Q = t N and c m are coefficients of the expansion in y. We search for a formal power series solution to this equation  , we obtain an overdetermined system of linear equations for the correlators C l 1 ...l n . From this system, one can actually pick certain equations in a clever way in specific order and obtain a recursive formula for the correlators. Suppose we wish to find the correlator C λ 1 ...λ •−1 λ • . Here λ • means the last part of the partition λ and we assume that indices of the correlator are sorted in descending order such that λ • is the smallest index (the correlator is of course symmetric as a function of its indices). We then consider to it (i.e., setting all times to zero after taking the partial derivative). We then get Here, μ denotes an ordered sequence of positive integer numbers (it can be an empty sequence), as opposed to a partition which is sorted sequence of positive integer numbers and l( μ) is the length of the sequence μ. # λ 1 is the number of parts of λ equal to 1, λ\λ • is partition λ without part λ • and λ\{λ • , ν} is partition λ with part λ • and all parts of ν removed.
The formula (5.4) is indeed a recursion: Correlators on all lines of the right hand side except the first one have sum of indices equal to |λ| − 2; the correlators on the first line, though having sum of indices equal to |λ|, all have a minimal index that is strictly less than λ • (thanks to the condition l( μ) ≥ 2).
Supplemented with the initial conditions C ∅ = 1 and C 1 = 0, this recursion allows us to calculate any particular correlator in a finite number of steps.
Example: The first step of the recursion looks like In the case of non-(q, t)-deformed Hermitean matrix model, the Virasoro constraints can be used to derive the so-calledŴ -operator representation of the partition function [24]. This derivation can be easily generalized to the case of Gaussian βensemble. However, at the moment we do not know how to generalize it to the (q, t)-Gaussian matrix model, as well as more complicated (q, t)-models. Existence ofŴ -operator form (the cut-and-join equation) is intimately linked with the existence of the spectral curve topological recursion for the model [7][8][9][10] and often signifies connections with enumerative geometry [25]. Pursuing these directions is definitely worth it in the future.

Gauge theory applications
From the point of view of gauge theories, the normalized correlators C l 1 ...l n are not very interesting. Instead, of primary interest are averages of Schur functions, which on the gauge theory side correspond to supersymmetric Wilson loops (see [12] and references therein). In the gauge theory on D 2 × S 1 , the expectation value of the Wilson loop along S 1 in representation λ corresponds to the following expression where s λ is Schur polynomial for the representation λ. If we choose N f = 2 with two opposite anti-fundamental masses ±m having the special values ±m f = ±q(1 − q) 1 2 and without fundamental contributions we can use the previously derived results. This choice of masses is by no means distinguished from the gauge theory point of view. Moreover, one could easily insert arbitrary values of the two anti-fundamental masses on the matrix model side as well, and this will not complicate q-Virasoro constraints at all (the factor 1 − q 2 (1 − q)y 2 in the first term of Eq. (5.1) will be replaced by (1 − m 1 y) (1 − m 2 y) but it will still remain quadratic, so non-Gaussian issues will still not arise). In what follows we present formulas at special values of the masses for conciseness.
Reexpanding the partition function in the basis of Schur functions (using Cauchy identity) it is straightforward to obtain constraints on Schur averages. In the case of the (q, t)-Gaussian model, the first few of them are: These expressions are examples of the expectation values of Wilson loops in D 2 × S 1 with the specific choices for the matter sector. In principle one may obtain the formulas for general values of fundamental mass by using the shifts in times. At the moment we do not know of a more direct way to write down the expressions of the form (6.3). Moreover, in the case of usual, non-(q, t)-deformed matrix models the Virasoro constraints look rather cumbersome in the Schur basis [26] so there is no reason to expect them to be simpler for the (q, t)-deformation. It is curious that from the point of view of (q, t)-deformed matrix models it is natural to study averages not of Schur functions, but rather of Macdonald functions M λ ({ p k }) [6]. Computer experiments show that in many examples Macdonald averages have a very simple, explicit form. For instance, for (q, t)-Gaussian matrix model one has though at the moment we do not know how to prove this formula. But, at least in principle, Schur averages could be obtained from the vector of Macdonald averages by applying the inverse matrix of (q, t)-Kostka coefficients [27]. The techniques suggested in this work can be combined with the free field realization of the gauge theory on squashed S 3 (ω 1 ,ω 2 ) [12]. There exists the generating function for the expectation values of supersymmetric Wilson loops here we use the notations from [12]. Upon certain assumptions this generating function satisfies two commuting sets of q-Virasoro constraints. It is straightforward to repeat the analysis from Sect. 4.1 keeping in mind that there is a change of variables x = e 2πi X/ω 1 (or x = e 2πi X/ω 2 for another copy of q-Virasoro). We leave this analysis for future studies.

Conclusion
In this paper we showed how to derive the set of Ward identities which a (q, t)eigenvalue model satisfies, starting from an integral representation of it. The derivation is based on a clever insertion of a full q-derivative under the integral. Under certain mild assumptions (which are to be checked in each particular case), these Ward identities are nothing but q-Virasoro constraints, which before could be derived only if the free field realization (in particular, the corresponding screening current) for the model was known.
The advantage of the method we present here is that it can be applied to any (q, t)eigenvalue model directly -there is no need to know the free field representation beforehand. Hence, the situation with (q, t)-eigenvalue models now becomes very similar to the situation with β-ensembles, where Virasoro constraints can be derived both from free field realization and from insertion of full derivatives under the integral. We think that having these two complementary perspectives is useful, as each of them highlights different aspects of the eigenvalue models.
It is worth noting that from the q-derivative insertion method one quite naturally obtains more constraints than from free field realization. Namely, there is an additional constraint T −1 Z = 0. During its derivation, one has to demand that the coefficient in front of trX −1 vanishes, so that one stays in the usual setting of matrix models. On the free field representation side, this extra equation is satisfied if one restricts attention to highest weights with zero momentum and it is very easy to overlook this special case.
There are several directions of further research that we think are interesting: • Application to concrete matrix models: The q-Virasoro constraints (4.20) can be viewed as a definition of the matrix model. A particular integral representation is then just one instance (or "phase") that satisfies this definition. It would be very interesting to study alternative phases of several objects which are of crucial importance in mathematical physics and are known to satisfy q-Virasoro constraints, such as ABJ(M) model and Nekrasov functions. •Ŵ -operator representation: In the case of usual Hermitean Gaussian matrix model, theŴ -operator representation is a consequence of Virasoro constraints. We hope that the explicit form of q-Virasoro constraints (5.4) will be useful in derivingŴ -operator representation for the (q, t)-Gaussian model. TheŴ -operator representation would then readily imply connection to enumerative geometry. • Large N limit and related questions: For non-(q, t)-deformed matrix models, the Virasoro constraints can be used to derive the so-called loop equations, which then may be solved in the large N limit [28]. A particular solution, satisfying certain mild assumptions, takes the form of a certain universal recursion procedure, the spectral curve topological recursion. Information about the concrete matrix model is translated into the initial algebrogeometric data for the recursion-the spectral curve. The spectral curve topological recursion is related to the Givental formalism for Cohomological field theories. It would be very interesting to see what the (q, t)-spectral curve topological recursion and (q, t)-Givental formalism look like. • Relation to Macdonald operators: The q-difference operators that we insert under the integral look very similar to the celebrated Macdonald operators [13,Equation (31)]. By understanding this connection better one could, perhaps, better understand the integrability of (q, t)-eigenvalue models. We hope that these questions would inspire some great future research.

A Special functions
A summary of the special functions used throughout the paper. The definition of the q-Pochhammer symbol is valid for |q| < 1, and the finite q-Pochhammer is given by The θ function is defined by The where # λ j is the number of parts j in partition λ and l(λ) is the length of partition λ.
The q-integral is defined as vanish at the boundary the q-shift can be rescaled away (A.6) One can also produce the q-integral by insertion of the operator O ν (x i ), which has the desired pole structure provided the choice of contours in Fig. 1. This can be seen using (A.5): The q-exponent has (at least) two definitions (A.10)

B.1 First starting relation
The starting relation (4.6) can be checked as follows. For the left hand side Whereas the right hand side becomes thus agreeing with the left hand side.

B.2 Second starting relation
The second starting relation can be verified since the left hand side matches the right hand side (B.7)

C Verification of algebraic identity
For the algebraic identity the proof of this is as follows. Starting by considering the contour integral of the right hand side where the contour is a big circle with radius tending to infinity, we see that this has apparent singularities at z = ∞ and z = 1 x i . But taking the limit z → ∞ (or equivalently α = 1 z → 0) we see that and thus there is no contribution from the singularity at z = ∞. Then looking at the other residue at z = 1 x i we see that which matches the product on the left hand side. To reproduce this singularity structure of the right hand side, we need to take this term and multiply by N i=1 generating the algebraic identity.

D.1 First constraint
we then change the order of the x i -integration and the ω-integration on the right hand side so that the poles at ω = q/x i are not relevant. Instead the residues at ω = z, y i , 0 and ∞ contribute with opposite sign. (The contour integral over all poles can be moved since we are on a sphere and thus the result is zero. Therefore, if we trade the residues at ω = q/x i the residues at ω = z, y i , 0 and ∞ will contribute with a minus sign.) This is illustrated in Fig. 3. Therefore Res ( Now the residues at ω = y i can be expanded in formal series since where in the last line we have used the Miwa variable t k = (1−t k ) i y k i . So these singularities do not contribute to the residue because we are interested in coefficients of {t k }.
Similarly we expect the relation (D.1) to hold at each order in a expansion around z → 0. Thus, there is no residue contribution at ω = z either. Then for the residue at ω = ∞ we consider and then use the parametrization ω = 1 α with dω = − 1 α 2 dα and instead take the limit α → 0. Now the contour will be anticlockwise but from the picture above we see that −Res(ω = ∞) = − (+Res(α = 0)) thus contributing with the same sign. Then we have Here we assumed c q (α) c q (qα) α=0 ∼ α 0 . This assumption is reasonable because if the ratio produced positive powers of α then the integrals would be vanishing, and if there were negative powers of α then Eq. (D.5) would get powers in z which would destroy the q-Virasoro constraints. Thus, we can impose that for the q-Virasoro constraints to hold we must have Then the equation becomes (D.6)

D.2 Second constraint
We have (D.7) Changing the order of integration similar to the first constraint equation , and considering the residue at ω = ∞ of the right hand side using the parametrization ω = 1 α with dω = − 1 α 2 dα and instead taking the limit α → 0, again assuming c q (α) c q (qα) α=0 ∼ α 0 . Using the above the equation becomes (D.9)

E Action of T 0 and T −1 from free field representation
The eigenvalue equations for T 0 and T −1 can also be seen from the free field representation. Starting with T 0 we use the explicit form of the generators T n of [12] (reproduced in (3.6)) using which we have where in the second line the sum over k drops out because all the positive modes in A k>0 annihilate the vacuum, and in the last line we choose the subset of zero momentum states with α = 0. For T −1 , we again use the explicit form of the generators using which we have where in the last line we set again α = 0. Thus, we recover the two eigenvalue equations from the free field representation in the case of zero momentum.