How can we prove that the expression
$$\gamma_{m k}(\lambda)=\sum_{i=1}^{m}\left(\lambda_{m i}+m-1\right)^{k} \prod_{j \neq i}\left(1-\frac{1}{\lambda_{m i}-\lambda_{m j}}\right)$$
is a polynomial in the variables $\lambda_{m1},...,\lambda_{mm}$, for all $m \geq 1 $ and $ m\geq k \geq 1$?
For a bit of context, this expression appears when studying Gelfand-Tsetlin representations in representation theory of Lie algebras, and it stands for the eigenvalues of operators $c_{mk}$ (which are the generators of the Gelfand-Tsetlin subalgebra of the enveloping algebra $\mathfrak{U}(\mathfrak{gl}(n))$) acting on a $\mathfrak{gl}(n)$-module $L(\lambda)$, when given in terms of some integer entries $\lambda_{mk}$. But it appears it can be proven without heavy machinery, although I tried and had no success.
Any ideas?
For a fixed $\,m=n\,$, let $\,\lambda_i = \lambda_{ni}\,$ and let $\displaystyle\,\sigma_k(\lambda) = \sum_{i=1}^{n} \lambda_{i}^{k} \prod_{j \neq i}\left(1-\frac{1}{\lambda_{i}-\lambda_{j}}\right)\,$. Since $\,\gamma_{mK}(\lambda)\,$ is a linear combination of $\,\sigma_k(\lambda)\,\big|_{\,k=0,1,\dots,K}\,$, it is enough to prove that $\,\sigma_k(\lambda)\,$ is a polynomial in $\,\lambda_i\,\big|_{\,i=1,2,\dots,n}\,$ for $\,k \ge 0\,$. Since $\,\sigma_k(\lambda)\,$ is symmetric in $\,\lambda_i\,$, this is equivalent to proving that it is a polynomial in the elementary symmetric polynomials i.e. the coefficients of $\displaystyle\,p(z)= \prod_{i=1}^n(z-\lambda_i)\,$.
Using that:
$$ \prod_{j \neq i}\left(1-\frac{1}{\lambda_{i}-\lambda_{j}}\right)=\frac{\prod_{j \neq i} (\lambda_{i}-1-\lambda_{j})}{\prod_{j \neq i} (\lambda_{i}-\lambda_{j})}=\frac{-\prod_{j} (\lambda_{i}-1-\lambda_{j})}{\prod_{j \neq i} (\lambda_{i}-\lambda_{j})}=- \frac{p(\lambda_i-1)}{p'(\lambda_i)} $$
it follows that:
$$ \sigma_k = -\sum_{i=1}^{n} \frac{\lambda_{i}^{k} \,p(\lambda_i-1)}{p'(\lambda_i)} $$
The sum reminisces of a Lagrange interpolation, and suggests looking at the polynomial:
$$ f(z) = \sum_{i=1}^{n} \frac{p(z)}{z-\lambda_i} \cdot\frac{\lambda_{i}^{k+1} \,p(\lambda_i-1)}{p'(\lambda_i)} $$
Since $\displaystyle\,\lim_{z \to \lambda_i} \frac{p(z)}{z-\lambda_i}=p'(\lambda_i)\,$ it follows that $\,f(\lambda_i) = \lambda_{i}^{k+1} \,p(\lambda_i-1)\,$, so $\,f(z) = z^{k+1}p(z-1)\,$ at the $\,n\,$ points $\,z=\lambda_1,\lambda_2,\dots, \lambda_n\,$.
Since $\,p(\lambda_i)=0\,$ that means $\,f(z)\,$ also coincides with $\,r_k(z) = z^{k+1} p(z-1) \bmod p(z)\,$, and must in fact be identical to it since $\,\deg r_k \lt \deg p = n\,$ and the two are equal at $\,n\,$ points.
Let $\,z^{k+1} p(z-1) = q_k(z) p(z) + r_k(z)\,$ where $\,q_k\,$ is a monic polynomial of degree $\,k+1\,$, then:
$$ z^{k+1} p(z-1) - q_k(z) p(z) \,=\, r_k(z) \,=\, f(z) \,=\, \sum_{i=1}^{n} \frac{p(z)}{z-\lambda_i} \cdot\frac{\lambda_{i}^{k+1} \,p(\lambda_i-1)}{p'(\lambda_i)} $$
$$ \implies \quad - q_k(0) p(0) \,=\, f(0) \,=\, -p(0)\,\sum_{i=1}^{n} \frac{\lambda_{i}^{k} \,p(\lambda_i-1)}{p'(\lambda_i)} \,=\, p(0)\,\sigma_k $$
$$ \implies\quad \sigma_k = -q_k(0) $$
Since $\,p(z)\,$ is monic, the coefficients of $\,q_k(z)\,$ are in the same ring with the coefficients of $\,p(z)\,$, see here for example. It then follows that $\,\sigma_k=-q_k(0)\,$ is a polynomial in the coefficients of $\,p(z)\,$.