Non-stationary difference equation and affine Laumon space III : Generalization to $\widehat{\mathfrak{gl}}_N$

In a series of papers we have considered a non-stationary difference equation which was originally discovered for the deformed Virasoro conformal block. The equation involves mass parameters and, when they are tuned appropriately, the equation is regarded as a quantum KZ equation for $U_q(A_{1}^{(1)})$. We introduce a $\widehat{\mathfrak{gl}}_N$ generalization of the non-stationary difference equation. The Hamiltonian is expressed in terms of $q$-commuting variables and allows both factorized forms and a normal ordered form. By specializing the mass parameters appropriately, the Hamiltonian can be identified with the $R$-matrix of the symmetric tensor representation of $U_q(A_{N-1}^{(1)})$, which in turn comes from the 3D (tetrahedron) $R$-matrix. We conjecture that the affine Laumon partition function of type $A_{N-1}^{(1)}$ gives a solution to our $\widehat{\mathfrak{gl}}_N$ non-stationary difference equation. As a check of our conjecture, we work out the four dimensional limit and find that the non-stationary difference equation reduces to the Fuji-Suzuki-Tsuda system.