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] 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 t1,t2 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 [2]. 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 [3]. 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 [4] and are naturally labelled by tuples of Young diagrams. The Hilbert space of the 2d theory is a representation of W N -algebra, 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 sec. 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 sec. 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 [5], obtain certain Baxter TQ equation. A particular case of this Baxter equation already appeared in the literature [6,7] as an equation determining the eigenstates of the integrals of motion associated to the Ding-Iohara-Miki algebra U t1,t2 ( 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 t1,t2 ( gl 1 ). It is natural to associate this representation with the states of a 3d field.
Finally, in sec. 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 Heisenberg 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) .
2 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 [8]). 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 data, since each solution is uniquely determined by the initial conditions 1 . In the mechanical systems the quantization amounts to replacing the coordinates on the space of solutions (e.g. initial coordinates and momenta) with noncommuting operators. In the field theory this procedure is the same in principle, though the space under consideration is infinite-dimensional, 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 .
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 they are killed by the boundary conditions on ∂D 2 . 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}: 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 (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: 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 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 11) 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: The corresponding function f (r) satisfies and thus 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 2 . 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. The coordinates on the space of solutions are the coefficients in the expansion: The Hilbert space H 3d tot is again 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: The pattern is clear, though explicit combinatorics of the mapping partitions to polynomials of a −(j,m) will not concern us here.
The 3d partition function (with m = 0 states removed) is thus equal to the product of two MacMahon functions 3 : (22) Adding back the solutions, which we have discraded corresponds to multiplying the partition function by an additional eta-function.
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 [10] (see also the review [9]). 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 4 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 [11]. 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. (23) 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. (23), but the argument is similar in both cases. Notice the difference between the indices in w n,m and W n,m : in the first case m ∈ Z ≥0 , 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 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: 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 (21) 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) : However, before applying Eq. (26) 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 [6,7]. Their solutions are labelled by 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" [12] obtained from refined topological strings [13]. 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. 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 5 : where is the matrix model measure. The difference operator acts on the measure as follows: so that half of the factors cancel with the extra factors in Eq. (29). The identity for the averages following from Eq. (29) We can use one more standard trick [14,15] 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 .
This equation can be viewed as a three-parametric (t 1 , t 2 and q) deformation of the Baxter equation associated to DIM algebra [7]. Notice also that the equation is symmetric under permutation of t 1 , t 2 and q t1t2 . 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 t1t2 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 (35) is then understood as an equation determining the saddle points x. The average signs in Eq. (35), therefore, disappear and it reduces to the DIM Baxter TQ equation: where Q(ξ) is a polynomial with roots x i . The only difference between our equation and the standard one considered in the literature [7] is that the roots of the polynomials K + (ξ) and K − (ξ) are a priori not related to each other. In [7], 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 (37) 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 [6], though the end results are different for the reason discussed above.
1. τ → −∞. Equations (38) simplify into: Suppose for a moment that K = 1, N = 2 and one of the variables, say x 1 sits at point z 1 . Then, equation (39) 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 We use a prime to distinguish partitions π ′(a) from the partitions π (a) from Eq. (28) 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.
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 z 1 (z 1 ) ⊗ · · · ⊗ M w K z K (z K ) with 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.
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 t1,t2 ( 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.

Triple Macdonald polynomials
In this section we build the states corresponding to the solutions of Bethe equations (39) 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 t1,t2 ( 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.
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 6 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 this space is written as follows: where C 0 is a small contour around zero and The equation for the eigenfunctions reads The eigenfunctions M (t1,t2) Y ( u|p a,n ) are called generalized Macdonald polynomials [14], [16] and depend on L-tuple of Young diagrams Y and L-tuple 7 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 nullvectors 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 [17]. 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 Rmatrix for the DIM algebra [18] 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 [18] for details): For x = t 1 t 2 the denominator in Eq. (46) 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 the R-matrix one can project out all the states |λ ⊗ |µ in the tensor product F (u) ⊗ F ( t1t2 x ), except those satisfying λ ⊃ µ. 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 (t1,t2) Y ( u|p a,n ) does not depend on L as long as L ≥ m. In particular, M (t1,t2) Y ( 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 (t1,t2) Y (u(t 1 t 2 ) 1−a |p a,n ). We denote the resulting polynomials by M (t1,t2,(t1t2) −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 (t1,t2,(t1t2) −1 ) π (p). The eigen-equation for triple Macdonald polynomials reads: Notice how the eigenvalues in the r.h.s. of Eq. (47) can be expressed as the sum of the Bethe roots x i from Eq. (42). 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: 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 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 (35) is what is called in [7] a module of "finite type", whereas the generator given in [20] 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) [13,21]. Moreover, the matrix model (27) 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 [22]. 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 [23]. It would be interesting to understand the relation between our results an those of [23].