3d field theory, plane partitions and triple Macdonald polynomials

We argue that MacMahon representation of Ding-Iohara-Miki (DIM) algebra spanned by plane partitions is closely related to the Hilbert space of a 3d field theory. Using affine matrix model we propose a generalization of Bethe equations associated to DIM algebra with solutions also labelled by plane partitions. In a certain limit we identify the eigenstates of the Bethe system as new triple Macdonald polynomials depending on an infinite number of families of time variables. We interpret these results as first hints of the existence of an integrable 3d field theory, in which DIM algebra plays the same role as affine algebras in 2d WZNW models.


Introduction
Nekrasov functions [1][2][3] obtained by applying the localization technique to supersymmetric gauge theories in four, five and six dimensions effectively sum up instanton contributions to their partition functions. Instanton configurations fixed under toric action of U(1) 2 t 1 ,t 2 are labelled by tuples of Young diagrams, which technically arise as ideals in the ring of polynomials in two variables, or as residues of the LMNS integral [4][5][6][7]. Young diagrams naturally remind one of the basis in the Hilbert space of a 2d field theory. It turns out that this connection can indeed be made precise leading to the famous AGT correspondence [8][9][10]. Each term in the instanton sum over partitions corresponds to a particular state in the Hilbert space (Verma module) of a 2d CFT. The states featuring in the correspondence are eigenstates of a certain family of CFT integrals of motion [11][12][13] and are naturally labelled by tuples of Young diagrams. The Hilbert space of the 2d theory is a representation of W Nalgebra, and is identified with the equivariant cohomology of the gauge theory instanton moduli space.
In this paper we develop the basics of a similar, though not yet as precise, identification for plane partitions, i.e. 3d Young diagrams. To this end in section 2 we study 3d free scalar field theory and identify (chiral half of) its Hilbert space as the space of plane partitions. We then consider in section 3, instead of the LMNS integral, the affine matrix model, which poles are labelled by plane partitions. We derive the quantum spectral curve for this model and, taking its Nekrasov-Shatashvili (NS) limit [14], obtain certain Baxter TQ JHEP06(2019)012 equation. A particular case of this Baxter equation already appeared in the literature 1 [16][17][18][19][20][21][22] as an equation determining the eigenstates of the integrals of motion associated to the Ding-Iohara-Miki algebra U t 1 ,t 2 ( gl 1 ). In this restricted case the solutions were labelled by K-tuples of Young diagrams, so that the eigenstates spanned a tensor product of K Fock representations (K 2d fields). However, the slight generalization of the Baxter equation obtained from the affine matrix model has solutions labelled by K-tuples of 3d Young diagrams. The states corresponding to these solutions, therefore, span the tensor product of what is called the MacMahon representations of U t 1 ,t 2 ( gl 1 ). It is natural to associate this representation with the states of a 3d field.
Finally, in section 4 we give a concrete construction of a basis in the MacMahon representation of DIM as a subrepresentation in an infinite tensor product of Fock representations. The basis, which we call the family of triple Macdonald polynomials, is formed by the eigenstates of the Cartan subalgebra of DIM algebra. This basis can also be obtained in a certain limit of DIM Baxter equation mentioned above. It is tempting to speculate that this construction implies the existence of an integrable 3d field theory, though we are currently very far from giving its concrete description. Ideologically, this connection might be justified by noticing that loop algebras (affine, Virasoro or W N ) naturally appear in 2d CFT, while DIM algebra is essentially a double loop algebra, and thus should be associated to a 3d theory. We use the standard notation for Young diagrams Y = [Y 1 , Y 2 , . . . , Y m ] with Y 1 ≥ Y 2 ≥ · · · ≥ Y m throughout the paper. Plane partition π = [π (1) , π (2) , . . . , π (n) ] is built from the "non-increasing" sequence of Young diagrams π (1) ⊃ π (2) ⊃ · · · ⊃ π (n) .

3d massless free scalar and 3d partitions
In this section we demonstrate how the Hilbert space of the simplest possible 3d field theory, the massless free scalar, is related to the space of plane partitions. We start by recalling the prototypical example of 2d free scalar, which Hilbert space factorizes into two chiral parts, each of which is a Fock space F with basis labelled by ordinary partitions (Young diagrams). We then analyze the 3d setup using similar arguments. In particular, we introduce an analogue of the holomorphic factorization and observe that the Hilbert space of the free scalar in 3d (after certain reduction) factorizes into a product of two MacMahon spaces spanned by plane partitions.
Our considerations are very elementary and the free scalar is only a toy model for more interesting 3d field theories (see e.g. the recent work [23]). Still, we believe that the general structure of the Hilbert space will remain the same as in the free case, especially if one considers integrable 3d models.

Warm-up: 2d free scalar
One way to quantize a classical model is to quantize the space of solutions to its classical equations of motion. One can also interpret it as the quantization of the space of initial JHEP06(2019)012 data, since each solution is uniquely determined by the initial conditions. 2 In the mechanical systems the quantization amounts to replacing the coordinates on the space of solutions (e.g. initial coordinates and momenta) with non-commuting operators. In the field theory this procedure is the same in principle, though the space under consideration is infinitedimensional, which may require more subtle treatment. However, for the free fields, no technical problems arise and the quantization can be carried out straightforwardly.
We will work in the Euclidean signature and the radial coordinate will play the role of time. The classical equations of motion for a 2d massless free scalar read Changing to the radial coordinates (r, ϕ) we get One can find a basis of solutions by separating the variables, i.e. φ(r, ϕ) = f (r)g(ϕ). From the fact that g(ϕ) needs to be single-valued we get The equation for f (r) is then and thus f (r) = r ±m . (2.5) At this point we need to specify on which part of the 2d plane the theory lives. We assume it to be a large disk D 2 around r = 0. Therefore, we can only consider positive powers of r, and thus f (r) = r |m| . We ignore the zero mode solutions f (r) = 1 and f (r) = ln r corresponding to m = 0. One can assume that the logarithmic mode is killed by the boundary conditions on ∂D 2 , while the constant mode could be discarded is one additionally requires f (r) to vanish at r = 0. 3 We will observe a similar situation when we come to the 3d case: part of the solutions will need to be discarded. One can also notice that the space of solutions which we have thrown away is exactly the space of solutions to the free scalar equations of motion in one dimension less. Indeed the functions 1 and ln r correspond to initial position and momentum of a particle on a line. In string theory this coordinate and momentum correspond to the motion of the string center of mass, while other solutions represent the harmonics. Regular solutions to the free field equations of motion are thus labelled by m ∈ Z\{0}:

JHEP06(2019)012
To quantize the space of solutions we need to introduce the coordinates. Those are the coefficients a −n andā −n in the expansion of a general field φ(r, ϕ) in the basis (2.6): They constitute in fact only half of the coordinates -the so-called positive-frequency parts. The other half corresponds to the solutions which are singular at r = 0, but are instead regular at r = ∞. If we introduce the complex structure on the 2d plane, so that the complex coordinates read z = re iϕ ,z = re −iϕ , then the coordinates a −n andā −n correspond to holomorphic and antiholomorphic fields respectively. The operators corresponding to the coordinates a −n andā −n commute between themselves and are the creation operators. They can be used to build all states of the model out of the vacuum state |∅ . After the choice of complex structure is made, the Hilbert space H 2d tot of the model, therefore, splits into a tensor product H 2d ⊗H 2d of holomorphic and antiholomorphic states: We can introduce the grading L 0 on H by assigning degree m to a −m . This choice of the degree corresponds to the action of the dilatation operator z∂ z on the solutions. Similar gradingL 0 z∂z acts onH 2d . The states of H 2d are labelled by Young diagrams {m 1 , m 2 , . . . , m k } (notice that the sequence is non-increasing). The partition function of the model can be written as a product of two Dedekind eta-functions: (2.10) where Y and W are 2d Young diagrams and |Y | denotes the total number of boxes in Y . The coordinates of the zero modes, that we have thrown away earlier (hence a prime in Z 2d ) would have contributed to the partition function by an overall factor of the form q P , which can be thought of as the partition function of a particle with a given momentum, or as the contribution of the center of mass coordinate. Let us recapitulate the result in the 2d case. The solutions of the equations of motion can be split into two parts: those regular at r = 0 and those regular at r = ∞. The first subspace corresponds to the creation operators. There is also a small subspace of zero modes, which needs to be discarded (or accounted for separately). After introduction JHEP06(2019)012 of the complex structure, creation operators further split into holomorphic and antiholomorphic subsets, each of which produces the Hilbert space with basis labelled by Young diagrams. The partition function is the product of two generating functions for the number of Young diagrams.

3d free scalar
We will now consider the 3d free scalar and analyze its quantization in the same spirit as the 2d case in the previous section. The equation of motion reads or in the spherical coordinates: where ∆ S 2 is the Laplacian operator on a unit sphere. Separating the variables we get φ(r, θ, ϕ) = f (r)g(θ, ϕ), where g(θ, ϕ) should be an eigenfunction of ∆ S 2 . It is thus given by the spherical harmonics: 14) The corresponding function f (r) satisfies r∂ r (r∂ r + 1)f (r) = j(j + 1)f (r), (2.15) and thus f (r) = r j or r −1−j (2.16) Only the first choice is regular at r = 0, so the solutions are given by Similarly to the 2d case we would like to discard certain solutions and split the remaining ones using some analogue of the complex structure. 4 The natural choice of the solutions to be discarded are those with m = 0. One way to justify this choice is to require the solutions to vanish at the t axis. This condition can be enforced on the distant boundary of the ball B 3 . There is only one discarded solution for any given j, so they can be thought of as corresponding to a 2d chiral scalar. Thus, the discarded solutions follow the pattern established in 2d: they correspond to the states of the field in one dimension less than the original model. We split the remaining solutions into those with positive and negative m calling them, by gross abuse of the terminology, holomorphic and antiholomorphic respectively. This 4 The relevant mathematical structure here is probably the contact structure represented by the standard one-form λ = dt + xdy − ydx = dt − i 2 (zdz − zdz) in R 3 . Notice that dλ = 2dx ∧ dy = idz ∧ dz represents the complex (or symplectic) structure on the constant t slices of R 3 .

JHEP06(2019)012
choice of splitting may be thought of as an additional geometric structure, which the model is supplemented with (see footnote 4). Just as in the 2d case, the equations of motion themselves do not imply any choice of splitting. However, in 2d an additional choice of complex structure allows one to define holomorphic and antiholomorphic modes. We would like to view the splitting of 3d modes in a similar fashion.
The coordinates on the space of solutions are the coefficients in the expansion: Let us mention that similar results for the quantization of the free scalar field in 3d can be obtained by the standard canonical quantization. To see this, we have to consider the Minkowski version of the free field model. More specifically, let the radial direction r become the time direction, while the space is spanned by (θ, ϕ), so that the sphere spacetime is locally R r × S 2 . In accordance with the general prescription we find the momentum corresponding to the field φ to be π = ∂ r φ. The canonical Poisson brackets are evaluated for a given time slice, e.g. r = 1: (2.20) are more conveniently written by expanding the field and its conjugate momentum in the complete basis of spherical harmonics Y j,m (θ, ϕ) on the sphere S 2 : In Euclidean field theory, the modes φ j,m , π j,m are recast into two sets of solutions to the field equations: those regular at r = 0 and those regular at r = ∞, as in eq. (2.16). In Minkowski field theory these two sets of modes are usually called the positive and negative frequency modes, because the powers of r become the frequencies after the Euclidean rotation ln r → i ln r. The coefficients of these modes become the classical versions of the creation and annihilation operators, which we denote by α † j,m and α j,m respectively:
The Hilbert space of the free field is obtained by acting with creation operators α † j,m on the vacuum vector |∅ . The splitting of the modes into "holomorphic" and "antiholomorphic" ones, which we have described above, corresponds to splitting the creation operators into those with m > 0 and those with m < 0 and throwing away those with m = 0. Thus we have the identification between the operators obtained from canonical quantization and The Hilbert space H 3d tot is (as in section 2.1) a tensor product of two factors H 3d ⊗H 3d . The space of holomorphic states H 3d is spanned by the vectors created from the vacuum |∅ by the strings of operators a −(j,m) : The grading on H 3d corresponds to the the total power of r in the solution, so that a −(j,m) has degree j.
Quite remarkably, the states on a given level in H 3d are in one to one correspondence with plane partitions:   From the table (2.32) one can see that at least on the first three levels the number of plane partitions and the number of states in H 3d coincide. The correspondence is in fact a combinatorial theorem 5 [26][27][28][29][30]. Let us briefly recall the steps of the proof: 1. Plane partitions are in one-to-one correspondence with pairs (P, Q) of Young tableaux sharing the same shape. Young tableaux of shape Y is a Young diagram Y filled with positive integers so that the numbers in the rows are non-increasing from left to right and the numbers in the columns are decreasing from top to bottom. An example of Young tableaux of shape [4, 2, 1] is Let π be a plane partition. It can be represented as a 2d array of non-negative integers π i,j , such that and π i,j = 0 for i or j sufficiently large. If we dissect π along the diagonal i = j, the cross section is a Young diagram Y defined by the non-increasing sequence of integers Y i = π i,i . The tableaux P and Q corresponding to π have the shape Y and are obtained from the "arm" and "leg" lengths corresponding to the boxes lying in the cross section, as shown in figure 1. The defining properties of the Young tableaux are satisfied for P and Q. Indeed, the numbers in the rows of P (resp. Q) are nonincreasing, since the corresponding arms (resp. legs) are stacked directly on top of each other, while the entries in a given column correspond to arms and lengths in a given horizontal slice of π, and therefore decrease as one goes further from the vertical coordinate axis.

JHEP06(2019)012
Young tableaux and generalized permutations, i.e. lexicographically ordered lists of pairs of integers. In this correspondence the first (resp. second) entry in each pair fills the tableaux Q (resp. P ). The detailed description of the correspondence and the proof that it is indeed one-to-one is given in [28,29] (see also [31]). In the example from figure 1 we get: 3. Any permutation σ of n objects can be written as an n × n matrix with entry σ ij equal to 1 if j = σ(i) and 0 otherwise. Similarly, any generalized permutation Σ can be written as an integer-valued matrix, in which Σ ij is given by the number of pairs i j in Σ. In our example, eq. (2.35), we get the matrix To an n×n matrix of integers Σ ij representing a generalized permutation we associate a state in H 3d : Notice that due to the properties of the RSK correspondence, the level of the state is equal to the total number of boxes in the initial partition, e.g. 7 in eq. (2.38) and figure 1.
Since the RSK correspondence is a bijection, so is our mapping of plane partitions to states in H 3d . The 3d partition function (with m = 0 states removed) is thus equal to the product of two MacMahon functions: 6 Adding back the solutions, which we have discraded corresponds to multiplying the partition function by an additional eta-function. The grading L 0 +L 0 can in fact be understood JHEP06(2019)012 as a scaling operator acting on the free boson field. The difference L 0 −L 0 doesn't seem to have such a direct physical interpretation.
What we see is that a particular splitting of the states of the free boson in 3d gives the analogue of holomorphic factorization for the 2d boson and that the states of the Hilbert space are labelled by (pairs of) plane partitions π.

MacMahon representation M c (u) of W 1+∞ and DIM algebras
The description of the chiral Hilbert space H 3d spanned by plane partitions π in terms of creation operators a −(j,m) might look a little unconventional. In this section we recall a very similar description of the MacMahon representation of the W 1+∞ algebra discovered in [25] (see also the review [24]). In this way we make a connection between the W 1+∞ algebra (and as we will see, DIM algebra in general) and the states of the 3d free scalar.
DIM algebra can be thought of as a t-deformation of the Lie algebra W 1+∞ . This Lie algebra is the central extension 7 of the algebra of difference operators W (n,m) q − nm 2 z n q mz∂z : where c 1 and c 2 are two central charges. One can equivalently understand W 1+∞ as the algebra of maps from the quantum torus (z, q z∂z ) to C. Notice that the relations of the algebra are manifestly invariant under the action of the SL(2, Z) automorphism group, which is the mapping class group of the (quantum) torus. In particular, the central charges transform as a doublet under this group. There is also a pair of grading operators (d 1 , d 2 ) which assign the weight (n, m) to the generator W n,m . We are not going to write down the relations for the DIM algebra for t = q, since they can be easily found in the literature [32,33]. There are, however, two general remarks. First of all, the t-deformed algebra is still SL(2, Z)-invariant. Secondly, the algebra is symmetric under the permutation of the parameters q, t −1 and t q . The remnant of this large symmetry is the symmetry of eq. (2.40) under q ↔ q −1 (one needs to simultaneously invert the sign of W n,m ).
We will describe the MacMahon representation M c (u) in the basis corresponding to the differential operators w n,m z n (z∂ z ) m , instead of the difference operators W n,m z n q mz∂z used in eq. (2.40), but the argument is similar in both cases. The difference between two choices of basis is only superficial. Indeed, one can expand the basis generators W n,m z n q mz∂z in ln q and get a series in the generators w n,m z n (z∂ z ) m . The two formulations of the algebra are completely equivalent, though the structure constants look much simpler in the basis W n,m . Notice the difference between the indices in w n,m and W n,m : in the first case m ∈ Z ≥0 since m labels the term in the Taylor expansion in ln q, while in the second n and m are general integers.
Consider the highest weight representation with central charge (c 1 , c 2 ) = (c, 0). The highest weight state |v , is the eigenstate of w 0,m annihilated by w n,m for n ≥ 1. The states JHEP06(2019)012 of the representation are obtained by acting on |v with the operators w n,m with n < 0 and m ≥ 0.
We further require there to be only a finite number of states on a given level. This is an extremely stringent requirement which fixes the allowed eigenvalues of w 0,m almost completely. To get a finite number of states on each level we need to have a lot of null-states. It turns out that the simplest nontrivial representation of this form has the null-subspaces generated by the following states [25]: where we have used the identification between the generators and the differential operators to write the null-states explicitly. Following the pattern of null states, we see that there are exactly n independent generators w −n,m for a given degree n > 0: those for which 0 ≤ m ≤ n − 1. All the other generators produce null-states. This situation is exactly parallel to the chiral free boson Hilbert space H 3d , which is generated by a −(j,m) with 1 ≤ m ≤ j. Thus, there is a direct equivalence between the holomorphic states of the 3d free boson (2.32) and the states of the MacMahon representation M c (u) of W 1+∞ : The only essential difference between the two pictures is that the creation and annihilation operators a −(j,m) with non-opposite vectors (j, m) commute while the commutator of two generators W n,m is nonzero in general.
Let us look for the analogy to this situation in the 2d free boson. Indeed, there are creation operators a −m , which commute for non-opposite m and there are also Virasoro generators L −n ∼ k∈Z : a k−n a −k : , which form a nontrivial Lie algebra. We can hypothesize that the same phenomenon happens in the 3d case and the W 1+∞ generators are expressed as bilinear combinations of bosonic operators a (j,m) : w n,m ? ∼ k,l : a n+k,m+l a −k,−l : . (2.43) However, before applying eq. (2.43) one needs to understand the limits of summation, the precise nature of the normal ordering and possible convergence issues. We will not attempt this task here.

Affine matrix models and DIM Bethe equations
In this section we derive the analogue of the quantum spectral curve for affine matrix models. Taking its Nekrasov-Shatashvili limit we obtain Baxter and Bethe equations, which are slight generalizations of those given in [16][17][18][19][20][21][22]. Their solutions are labelled by JHEP06(2019)012 plane partitions. Following the results of the previous section this hints at a possible connection with an integrable 3d field theory. Affine matrix model is one of the species from the zoo of "network matrix models" [34][35][36][37] obtained from refined topological strings [40][41][42][43][44][45]. Concretely, this model corresponds to the compactified strip geometry. In the geometric engineering terms it gives the 5d N = 1 U(K) gauge theory with an extra adjoint multiplet, or equivalently 6d abelian linear quiver. The matrix model average of a function f ( x) is given by: where Z is the integral without the insertion and the contour C can be chosen to encircle the poles of the integrand. The parameters z i and w j can be viewed as Kähler parameters of the strip geometry, or Coulomb moduli of the 5d theory. These poles correspond to K-tuples of plane partitions π (a) , a = 1, . . . , K: The partitions π (a) have the total floor area N , while the total number of boxes can be arbitrary. The residues at the poles transform naturally under arbitrary permutations of t 1 , t 2 and q −1 and the simultaneous transpositions of the plane partitions π (a) (assuming N to be sufficiently large).

Quantum spectral curve
We follow the standard technique to obtain the quantum spectral curve (or loop equations, or Ward identities, or qq-characters) for the model and consider an integral of a total difference: 8 is the matrix model measure. The difference operator acts on the measure as follows: The contour integral in the definition of the model is equivalent to a Jackson integral, for which total differences play the role of total derivatives.

JHEP06(2019)012
so that half of the factors cancel with the extra factors in eq. (3.3). The identity for the averages following from eq.
We can use one more standard trick [46,47] and rewrite the sums under the average as contour integrals over an auxiliary parameter y: where the contour C x encircles all the points x j . Deforming the contour we pick up the residues at y = ξ and at y = 0, ∞. The former residue gives the expression for the quantum spectral curve while the latter two produce polynomials in ξ with coefficients polynomially depending on x i . where (3.10)

JHEP06(2019)012
This equation can be viewed as a three-parametric (t 1 , t 2 and q) deformation of the Baxter equation associated to DIM algebra [20][21][22]. Notice also that the equation is symmetric under permutation of t 1 , t 2 and q t 1 t 2 . This symmetry combines with the symmetry of the integrand under the permutation of t 1 , t 2 and q −1 , which we mentioned in the beginning of this section. Overall, the matrix model is symmetric under permutation of four parameters, t 1 , t 2 , q −1 and q t 1 t 2 with their product being equal to one.

NS limit and Bethe equations
In the NS limit q → 1 the matrix model integral can be evaluated by saddle point method. The quantum spectral curve equation (3.9) is then understood as an equation determining the saddle points x. The average signs in eq. (3.9), therefore, disappear and it reduces to the DIM Baxter TQ equation: (3.11) where Q(ξ) is a polynomial with roots x i . The only difference between our equation and the standard one considered in the literature [20][21][22] is that the roots of the polynomials K + (ξ) and K − (ξ) are a priori not related to each other. In [20][21][22], however, one had the condition K − (ξ) = K + (t 1 t 2 ξ). In the matrix model language this constraint relating w a and z a is interpreted as the degeneracy condition for the vertex operators sitting at points z a -their weights (dimensions) take special, quantized values. As we will see momentarily, lifting of the constraint has dramatic consequences for the structure of the solutions: with the constraint they are enumerated (up to permutations) by K-tuples of Young diagrams, while without it the solutions are much more numerous and correspond to K-tuples of plane partitions.
Bethe equations following from the Baxter equation (3.11) are: where e τ = lim q→1 q u+1 (we assume that u scales as 1 ln q in the NS limit). Let us consider two limits τ → ±∞ of these equations and see how plane partitions arise as solutions. Our analysis is very similar to that of [16][17][18][19], though the end results are different for the reason discussed above.

JHEP06(2019)012
Suppose for a moment that K = 1, N = 2 and one of the variables, say x 1 sits at point z 1 . Then, equation (3.13) for i = 1 is already satisfied and we need not consider it anymore. The next variable, x 2 should solve one of the remaining equations. It cannot sit at z 1 because of the denominator (x 1 − x 2 ), and therefore should be either t 1 z 1 , t 2 z 1 or (t 1 t 2 ) −1 z 1 . These three solutions correspond to three plane partitions with N = 2 boxes.
In general one can see that the solutions are labelled by K-tuples of plane partitions π (a) with total of N boxes. They are given by (3.14) We use a prime to distinguish partitions π (a) from the partitions π (a) from eq. (3.2) labelling the poles of the affine matrix model integrand. In the NS limit the poles condense, so that the size of the saddle point partitions π (a) is actually infinite. Notice that N is the total number of boxes in π (a) , while the partitions π (a) have total floor area N and arbitrary number of boxes.
2. τ → +∞. This limit is analogous to the previous one: The solutions are again labelled by plane partitions, but the roots cluster around the points w a , not z a : The solutions for general τ extrapolate between two limits τ → ±∞, but are still enumerated by plane partitions. We have verified this claim by numerically solving the equations for small N and K. We can then ask what happens in the limit when w a → t 1 t 2 z a ? One can notice that in this case plane partitions with height more than one cease to be solutions, because of extra cancellations between the factors in K + and K − . We have also verified this effect numerically. More generally, if w a → (t 1 t 2 ) n z a the solutions are enumerated by plane partitions with height at most n.
The solutions to Bethe equations correspond to the basis of eigenstates of DIM integrals of motion in the corresponding DIM representation. For general K + and K − the representation is a tensor product of MacMahon modules M w 1 spectral parameters z a and central charges (ln wa za , 0). The basis is given by the K-tuples of plane partitions. For special degenerate value of the central charge equal to (ln(t 1 t 2 ), 0) the MacMahon module becomes reducible. After factoring out the invariant subspace one obtains the Fock representation with standard basis labelled by Young diagrams. The disappearing solutions of the Bethe equations correspond to the invariant subspaces factored out from the MacMahon modules. For w a → (t 1 t 2 ) n z a the MacMahon module degenerates into a subquotient of tensor product of n Fock modules. Correspondingly, the central charge of the resulting representation is (n ln(t 1 t 2 ), 0).

JHEP06(2019)012
Let us summarize the results of this section. We have derived the quantum spectral curve equation for the affine matrix model and studied the Bethe equations obtained in the NS limit. The solutions of these equations are in general labelled by K-tuples of plane partitions and span the tensor product of MacMahon representations of DIM algebra U t 1 ,t 2 ( gl 1 ). The Bethe equations, therefore describe certain (eigen)states of a 3d field field theory. In the next section we give an explicit construction of the basis of eigenstates in the single MacMahon module.
Let us emphasize that the Bethe eigenstates for a basis in a representation of DIM algebra, i.e. there is an action of DIM algebra on them. The action is described in terms of combinatorics of 3d Young diagrams, i.e. the currents generating the DIM algebra act on a Young diagram by either adding or removing a box with some weight depending on the diagram and the box. However, it seems that this action is not easy to see in the matrix model picture (3.1). The toy example of a similar situation is the XXX spin chain where Bethe eigenstates form a special basis in the representation C 2 ⊗ · · · ⊗ C 2 of the Yangian Y [sl 2 ], though the explicit action on them is somewhat complicated. It would be interesting to work out the DIM action on Bethe eigenstates using the Ward identities of the matrix model (3.1).

Triple Macdonald polynomials
In this section we build the states corresponding to the solutions of Bethe equations (3.13) for the simplest case of K = 1. We model the MacMahon module as a subspace inside the infinite tensor product of Fock modules F (u) of DIM algebra U t 1 ,t 2 ( gl 1 ) with specially adjusted spectral parameters. The infinite tensor product of Fock spaces, each of which is effectively a 2d free boson, can be seen as a discretized construction of a 3d field. One can view such a discretization as follows. The Fock space F (u) is the Hilbert space of a 2d boson on a circle S 1 (the time direction is assumed to be R). An infinite tensor product of Fock spaces is then the Hilbert space of an infinite family of free 2d bosons on a collection of circles. This infinite collection of independent bosons on circles can be seen as an infinite set of modes of a single 3d boson on a 2d space slice e.g. on a torus S 1 × S 1 . The modes in this case naturally have two discrete indices, corresponding to the momenta along two circles. As we have argued in section 2.3, in the case of MacMahon representation the relevant 2d space slice is S 2 rather than S 1 × S 1 . However, this only changes the mode labels, while the overall construction of "3d from 2d" still goes through, as we show below.
Concretely, we consider the eigenfunctions of the Cartan subalgebra of DIM algebra. The first element of this subalgebra is the zero mode of the current x + 0 , which corresponds to generator W 0,1 in the W 1+∞ notation. We will diagonalize this element in the tensor product of Fock representations 9 F (u 1 ) ⊗ · · · ⊗ F (u L ), which we identify with the space of polynomials in L families of time variables p a,n , a = 1, . . . , L, n ≥ 1. The action of x + 0 on 9 The definition of the tensor product of representations requires the choice of the coproduct ∆ in the DIM algebra. This amounts to choosing a preferred direction in the space of central charge vectors. Here we choose the horizontal direction.

JHEP06(2019)012
this space is written as follows: where C 0 is a small contour around zero and (4. 2) The equation for the eigenfunctions reads The eigenfunctions M ( u|p a,n ) are called generalized Macdonald polynomials [46,48,49] and depend on L-tuple of Young diagrams Y and L-tuple 10 of spectral parameters u.
For certain special values of spectral parameters called the resonances, i.e. when u a = u a+1 t m 1 t l 2 , part of the tensor product of Fock modules completely decouples (i.e. the corresponding states become null-vectors inside the DIM representation) leaving only polynomials corresponding to the special L-tuples of Young diagrams. The same phenomenon can be observed when considering the Nekrasov expansion of the conformal block. To get the right answer for degenerate values of the intermediate dimensions one should constraint the Young diagrams appearing in the expansion to satisfy the Burge conditions [50][51][52]. We would like to consider here the simplest case of the resonance when u a = u(t 1 t 2 ) 1−a . In this case the surviving L-tuples of Young diagrams form a plane partition π = {Y 1 , . . . , Y L } of height L.
Another way to see the decoupling of states for spectral parameters in resonance is to consider the R-matrix for the DIM algebra [53][54][55][56] acting in the tensor product of Fock spaces F (u) ⊗ F ( u x ). The R-matrix in the resonant case becomes singular. Indeed, the eigenvalues of the R-matrix can be read off from its expression in the vertical representation (see [53][54][55][56] for details): .
(4.4) For x = t 1 t 2 the denominator in eq. (4.4) is finite for any Young diagrams λ and µ while the numerator vanishes whenever µ doesn't fit inside λ. Since the operator x + 0 is an element of DIM algebra, its action on the tensor product commutes with the R-matrix. Thus, using 10 More precisely M (t 1 ,t 2 ) Y ( u|pa,n) depend only on (L − 1) ratios of ua.
The effect here is similar to the fusion of spins in a spin chain. Take as an example the rational R-matrix R(x − y) = 1 − P x−y acting in the tensor product of two fundamental evaluation representations C n at points x and y (P is the permutation operator). This R-matrix becomes a projector for the resonant values of the spectral parameters x = y ± 1. For these values one can take a projection of the tensor product using R(±1) and, iterating this procedure, build evaluation representations with general spin.
Back to the DIM case: we notice that generalized Macdonald polynomials are stable in the following sense. Consider an L-tuple of Young diagrams Y , in which only first m diagrams are non-empty. Then M (t 1 ,t 2 ) Y ( u|p a,n ) does not depend on L as long as L ≥ m. In ( u|p a,n ) in this case depends only on the first m times p a,n , a = 1, . . . , m but not on those with m < a ≤ L. Using this stability property we can take the limit L → ∞ of the polynomials M (t 1 ,t 2 ) Y (u(t 1 t 2 ) 1−a |p a,n ). We denote the resulting polynomials by M (t 1 ,t 2 ,(t 1 t 2 ) −1 ) π (p) and call them triple Macdonald polynomials. They depend on the plane partition π and an infinite number of families of time variables p a,n a ≥ 1, though for any concrete π only a finite number of p a,n enters M (t 1 ,t 2 ,(t 1 t 2 ) −1 ) π (p). The eigen-equation for triple Macdonald polynomials reads: Notice how the eigenvalues in the r.h.s. of eq. (4.5) can be expressed as the sum of the Bethe roots x i from eq. (3.16). The eigenvalues are invariant with respect to the permutation of t 1 , t 2 and (t 1 t 2 ) −1 and the corresponding simultaneous transposition of π. Only part of this symmetry exchanging t 1 and t 2 is manifest in our description, and it would be very interesting to understand the other hidden part. The most elementary property of Macdonald polynomials is that for t 1 = t −1 2 they reduce to the product of Schur polynomials: s π (i) (p i,n ) (4.6) They also satisfy the inversion relation analogous to that of ordinary Macdonald polynomials: where π T denotes one of the transpositions of π, which acts on the slices of fixed height, and ). (4.8)

JHEP06(2019)012
and the corresponding solutions to Bethe equations. We give arguments that the natural framework to work with such representations is not the 2d fields, but 3d ones, where the (chiral) Hilbert space is spanned by plane partitions.
Let us give directions in which we would like to extend our results in the future.
Matrix models of the type we considered can be understood as Dotsenko-Fateev representations [57,58] of certain W -algebra conformal blocks. The quantum spectral curve (3.9) corresponds to a particular generator of the W -algebra commuting with a set of screening charges, whose correlator gives the affine matrix model. This W -algebra is associated to a circular quiver with one node [59][60][61]. The advantage of our quantum spectral curve is that it is a finite expression, whereas the generators given in [59][60][61] were infinite series. One can notice that the origin of this difference lies in the difference between classes of DIM representations, since the W -algebra provides the quantization of DIM representation ring. The representation corresponding to (3.9) is what is called in [20][21][22] a module of "finite type", whereas the generator given in [59][60][61] is associated to a Fock module, which is not of finite type. It should be possible to write down the relations of the W -algebra associated to DIM algebra explicitly in the basis corresponding to finite type modules.
Plane partitions naturally arise in the crystal melting models of (refined) topological vertex C λµν (q, t) [38][39][40][41][42][43][44][45]62]. Moreover, the matrix model (3.1) actually reduces to the refined crystal melting partition function in the NS limit, i.e. q → 1 (the remaining parameters t 1 , t 2 become q and t −1 ). The only missing ingredients are the non-empty asymptotics of the plane partitions featuring in the melting crystal model. We plan to investigate this connection in the future. We would also like to point out one mysterious phenomenon along this direction. Refined topological vertex C λµν (q, t) can be identified with the intertwiner of Fock representations of DIM algebra [63]. Simultaneously, the Baxter equations obtained from the corresponding matrix model in the NS limit are relations in the representation ring of DIM algebra. How to make a direct connection between an intertwiner acting between states of the representations and a relation between the products of representation spaces?
Very recently a connection between solid (i.e. 4d) partitions and gauge theory was put forward in [64]. It would be interesting to understand the relation between our results an those of [64].