  Based on the framework of Plamen Iliev, multivariate Meixner polynomials are
constructed explicitly as Birth and Death polynomials. 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)\in\mathbb{N}_0^n$ are
$B_j(x)=\bigl(\beta+\sum_{i=1}^nx_j\bigr)$ and $D_j(x)=c_j^{-1}x_j$, $0<c_j$,
$j=1,\ldots,n$, $\sum_{j=1}^nc_j<1$. The corresponding stationary distribution
is
$(\beta)_{\sum_{j=1}^nc_j}\prod_{j=1}^n(c_j^{x_j}/x_j!)(1-\sum_{j=1}^nc_j)^\beta$,
the trivial $n$-variable generalisation of the orthogonality weight of the
single variable Meixner polynomials. The polynomials, depending on $n+1$
parameters ($\{c_i\}$ and $\beta$), 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 Meixner
polynomials. The polynomials are truncated $(n+1,2n+2)$ hypergeometric
functions of Aomoto-Gelfand. The polynomials and the derivation are very
similar to those of the multivariate Krawtchouk polynomials reported recently.
