Multivariate Kawtchouk polynomials as Birth and Death polynomials

Multivariate Krawtchouk polynomials are constructed explicitly as Birth and Death polynomials, which have the nearest neighbour interactions. They form the complete set of eigenpolynomials of a birth and death process with the birth and death rates at population $x=(x_1,\ldots,x_n)$ are $B_j(x)=\bigl(N-\sum_{i=1}^nx_i\bigr)$ and $D_j(x)=p_i^{-1}x_j$, $0<p_j$, $j=1,\ldots,n$. The corresponding stationary distribution is the multinomial distribution with the probabilities $\{\eta_i\}$, $\eta_i= p_i/(1+\sum_{j=1}^np_j)$. The polynomials, depending on $n+1$ parameters ($\{p_i\}$ and $N$), satisfy the difference equation with the coefficients $B_j(x)$ and $D_j(x)$ $j=1,\ldots,n$, which is the straightforward generalisation of the difference equation governing the single variable Krawtchouk polynomials. The polynomials are truncated $(n+1,2n+2)$ hypergeometric functions of Aomoto-Gelfand. The divariate Rahman polynomials are identified as the dual polynomials with a special parametrisation.