Multivariate Meixner polynomials as Birth and Death polynomials

Ryu Sasaki

Department of Physics and Astronomy, Tokyo University of Science, Noda
278-8510, Japan

###### Abstract

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 $\boldsymbol{x}=(x_{1},\ldots,x_{n})\in\mathbb{N}_{0}^{n}$
are $B_{j}(\boldsymbol{x})=\bigl{(}\beta+\sum_{i=1}^{n}x_{j}\bigr{)}$ and
$D_{j}(\boldsymbol{x})=c_{j}^{-1}x_{j}$, $0<c_{j}$, $j=1,\ldots,n$,
$\sum_{j=1}^{n}c_{j}<1$. The corresponding stationary distribution is
$(\beta)_{\sum_{j=1}^{n}c_{j}}\prod_{j=1}^{n}(c_{j}^{x_{j}}/x_{j}!)(1-\sum_{j=1}^{n}c_{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}(\boldsymbol{x})$ and $D_{j}(\boldsymbol{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.

##  1 Introduction

As the second member of the discrete multivariate hypergeometric orthogonal
polynomials of Askey scheme [2, 20, 23, 30], i.e. those satisfying second
order difference equations with the nearest neighbour interactions, the
multivariate Meixner polynomials are constructed explicitly as the
multivariate Birth and Death (BD) [5, 21, 20] polynomials. Multivariate BD
processes are the nearest neighbour interactions of multi-dimensional discrete
systems, the most basic type of interactions, like the well known Ising
models.

The main strategy is essentially the same as that of the multivariate
Krawtchouk orthogonal polynomials reported recently [30]. The basic framework,
the orthogonality measure, the design of the truncated hypergeometric
functions and the generating function connecting them, is provided by a lucent
paper of Iiiev [18] which follows the trends of many preceding works, e.g. [8,
25]. The additional input is the reformulation of the birth and death problems
with the link to the associated self-adjoint matrices [27, 29, 28]. It is
markedly different from the traditional Karlin-McGregor approach [21, 20].
Since the central motivation is the pursuit of multivariate orthogonal
polynomials obeying second order difference equations, this paper is rather
distant from other works on multivariate orthogonal polynomials [3, 31, 33,
24, 26, 15, 16, 11, 14, 22, 12, 13, 17, 7, 32, 4, 10].

This paper is organised as follows. Starting with a short résumé of the single
variable Meixner polynomials in section 2, the necessary ingredients for the
construction of multivariate Meixner polynomials are listed in the logical
order in section 3. In section 4, Iliev’s framework for the multivariate
Meixner polynomials is reproduced in my notation. The main contents are in
section 5. An overview of multivariate birth and death problem is
recapitulated in section 5.1. The birth and death rates are introduced in
section 5.2. The orthogonality of the Meixner polynomials is proved in section
5.3. The second order difference equation of the multivariate Meixner
polynomials is demonstrated by using the generating function in section 5.4.
The solution of the BD problem is presented in section 5.5. The exceptional
cases are remarked in section 5.6.

##  2 Single variable Meixner polynomials

Let us start with the summary of the single variable Meixner polynomials [23]
with real positive parameters $1>c>0$, $\beta>0$,

$\qquad
P_{m}(\beta,c;x)\stackrel{{\scriptstyle\text{def}}}{{=}}{}_{2}F_{1}\Bigl{(}\genfrac{}{}{0.0pt}{}{-m,\,-x}{\beta}\Bigm{|}1-c^{-1}\Bigr{)}=\sum_{k\in{\mathbb{N}}_{0}}\frac{(-m)_{k}(-x)_{k}}{(\beta)_{k}}\frac{(1-\frac{1}{c})^{k}}{k!},\quad
m\in\mathbb{N}_{0},$ (2.1)

which are obtained by the generating function $G(\beta,c,x;t)$,

$\displaystyle
G(\beta,c,x;t)\stackrel{{\scriptstyle\text{def}}}{{=}}\left(1-\frac{t}{c}\right)^{x}(1-t)^{-\beta-x}=\sum_{m\in{\mathbb{N}}_{0}}\frac{(\beta)_{m}}{m!}P_{m}(\beta,c;x)t^{n},$
(2.2)

in which $(a)_{n}$ is the shifted factorial defined for $a\in\mathbb{C}$ and
nonnegative integer $n$, $(a)_{0}=1$, $(a)_{n}=\prod_{k=0}^{n-1}(a+k)$, $n\geq
1$. They satisfy the orthogonality relations with the normalised orthogonality
weight $W(\beta,c;x)$,

$\displaystyle\qquad
W(\beta,c;x)\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{(\beta)_{x}c^{x}}{x!}(1-c)^{\beta},$
(2.3)
$\displaystyle\sum_{x\in{\mathbb{N}}_{0}}W(\beta,c;x)P_{m}(\beta,c;x)P_{m^{\prime}}(\beta,c;x)=\frac{1}{(\beta)_{n}\frac{c^{m}}{m!}}\,\delta_{m\,m^{\prime}},$
(2.4)

and the second order difference equations are

$\displaystyle\widetilde{\mathcal{H}}\stackrel{{\scriptstyle\text{def}}}{{=}}(\beta+x)(1-e^{\partial})+\frac{x}{c}(1-e^{-\partial}),\qquad\partial=\frac{d}{dx},\quad
e^{\pm\partial}f(x)=f(x\pm 1),$ (2.5)
$\displaystyle\widetilde{\mathcal{H}}P_{m}(\beta,c;x)=\frac{1-c}{c}mP_{m}(\beta,c;x),\qquad
m\in{\mathbb{N}}_{0}.$ (2.6)

The eigenvalue is easily guessed as

$\widetilde{\mathcal{H}}x^{m}=m\left(-1+\frac{1}{c}\right)x^{m}+\text{lower
degrees},\qquad m\in{\mathbb{N}}_{0}.$

It is straightforward to verify (2.6) in terms of the generating function,

$\widetilde{\mathcal{H}}G(\beta,c,x;t)=\frac{1-c}{c}t\frac{\partial}{\partial
t}G(\beta,c,x;t).$ (2.7)

##  3 Path to multivariate Meixner polynomials

In this paper I present the explicit forms of multivariate Meixner polynomials
satisfying second order difference equations in a similar way to the single
variable Meixner polynomials shown above. The necessary ingredients are

  1. 1.

(normalised) orthogonality weight $W(\boldsymbol{x})$

  2. 2.

generating function $G(\boldsymbol{x})$

  3. 3.

general form of the polynomial $P_{\boldsymbol{m}}(\boldsymbol{x})$

  4. 4.

proof of the orthogonality
$\sum_{\boldsymbol{x}\in{\mathbb{N}}_{0}^{n}}W(\boldsymbol{x})P_{\boldsymbol{m}}(\boldsymbol{x})P_{{\boldsymbol{m}}^{\prime}}(\boldsymbol{x})=0$,
$\boldsymbol{m}\neq\boldsymbol{m}^{\prime}$

  5. 5.

second order difference operator $\widetilde{\mathcal{H}}$

  6. 6.

eigenvalue spectrum $\mathcal{E}(\boldsymbol{m})$

  7. 7.

proof of
$\widetilde{\mathcal{H}}P_{\boldsymbol{m}}(\boldsymbol{x})=\mathcal{E}(\boldsymbol{m})P_{\boldsymbol{m}}(\boldsymbol{x})$

Plamen Iliev provided the first four ingredients in [18], which will be cited
as I. The formulation of multivariate Birth and Death polynomials introduced
by myself [30] provides the rest and the explicit forms of the multivariate
Meixner polynomials are obtained.

##  4 Framework for the multivariate Meixner polynomials due to Iliev [18]

###### Definition 4.1

The normalised orthogonality weight of the multivariate Meixner polynomials is
a simple $n$-variable generalisation of the single variable one (2.3)

$W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{(\beta)_{|x|}\boldsymbol{c}^{\boldsymbol{x}}}{\boldsymbol{x}!}(1-|c|)^{\beta},\quad\beta>0,$
(4.1)

in which the probability parameters $\\{c_{i}>0\\}$ are restricted by the
summability of $W$,

$\displaystyle{\boldsymbol{x}}$
$\displaystyle=(x_{1},x_{2},\ldots,x_{n})\in{\mathbb{N}}_{0}^{n},\quad|x|\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{i=1}^{n}x_{i},\quad{\boldsymbol{x}}!\stackrel{{\scriptstyle\text{def}}}{{=}}\prod_{i=1}^{n}x_{i}!,$
$\displaystyle{\boldsymbol{c}}$
$\displaystyle=(c_{1},c_{2},\ldots,c_{n})\in{\mathbb{R}}_{>0}^{n},\,\quad|c|\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{i=1}^{n}c_{i},\quad{\boldsymbol{c}}^{\boldsymbol{x}}\stackrel{{\scriptstyle\text{def}}}{{=}}\prod_{i=1}^{n}c_{i}^{x_{i}},$
$\displaystyle\sum_{\boldsymbol{x}\in{\mathbb{N}}_{0}^{n}}W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})=1\quad\Rightarrow
0<|c|<1.$ (4.2)

First, prepare $n^{2}$ real parameters $\boldsymbol{u}=\\{u_{i\,j}\\}$,
$i,j=1,\ldots,n$ constrained by the conditions

$\displaystyle\sum_{i=1}^{n}c_{i}u_{i\,j}=|c|-1,\quad j=1,\ldots,n,$ (4.3)
$\displaystyle\sum_{i=1}^{n}c_{i}u_{i\,j}u_{i\,k}=|c|-1,\quad j\neq k,\quad
j,k=1,\ldots,n.$ (4.4)

By introducing another closely related $n^{2}$ real parameters
$\boldsymbol{b}=\\{b_{i\,j}\\}$,

$b_{i\,j}\stackrel{{\scriptstyle\text{def}}}{{=}}1-u_{i\,j},\qquad
i,j=1,\ldots,n,$ (4.5)

which are Iliev’s original parameters, these conditions look simpler

$\displaystyle\sum_{i=1}^{n}c_{i}b_{i\,j}=1,\quad j=1,\ldots,n,$ (I.2.3a)
$\displaystyle\sum_{i=1}^{n}c_{i}b_{i\,j}b_{i\,k}=1,\quad j\neq k,\quad
j,k=1,\ldots,n.$ (I.2.3b)

Similar conditions have been introduced by many authors for the construction
of multivariate orthogonal polynomials [8, 9, 25, 10, 30]. It should be
stressed that these $n(n+1)/2$ conditions are not enough to determine $n^{2}$
parameters ${\boldsymbol{u}}$ completely.

###### Definition 4.2

Iliev [I.(2.3b)] also introduces another set of probability parameters
$\\{\bar{c}_{j}\\}$, by

$1-|c|+\sum_{i=1}^{n}c_{i}u_{i\,j}^{2}=\frac{1-|c|}{\bar{c}_{j}},\quad
j=1,\ldots,n.$ (4.6)

###### Definition 4.3

The generating function $G$ is defined by (I.(2.6))

$G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})\stackrel{{\scriptstyle\text{def}}}{{=}}(1-|t|)^{-\beta-|x|}\prod_{i=1}^{n}\left(1-\sum_{j=1}^{n}b_{i\,j}t_{j}\right)^{x_{i}},$
(4.7)

in which

${\boldsymbol{t}}=(t_{1},t_{2},\ldots,t_{n})\in{\mathbb{C}}^{n},\quad|t|\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{i=1}^{n}t_{i}.$

###### Definition 4.4

The polynomials $P_{\boldsymbol{m}}(\boldsymbol{x})$ are defined by the
expansion of $G$ (I.(2.7)) around ${\boldsymbol{t}}={\boldsymbol{0}}$,

$G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})=\sum_{{\boldsymbol{m}}\in{\mathbb{N}}_{0}^{n}}\frac{(\beta)_{|m|}}{{\boldsymbol{m}}!}P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};{\boldsymbol{x}}){\boldsymbol{t}}^{\boldsymbol{m}},\quad{\boldsymbol{m}}=(m_{1},\ldots,m_{n})\in{\mathbb{N}}_{0}^{n},\quad{\boldsymbol{t}}^{\boldsymbol{m}}=\prod_{i=1}^{n}t_{i}^{m_{i}}.$
(4.8)

His two main Theorems are on the orthogonality and the concrete form of
$P_{\boldsymbol{m}}({\boldsymbol{x}})$.

###### Theorem 4.5

Orthogonality relation reads

$\displaystyle\sum_{{\boldsymbol{x}}\in{\mathbb{N}}_{0}^{n}}W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};{\boldsymbol{x}})P_{{\boldsymbol{m}}^{\prime}}(\beta,{\boldsymbol{u}};{\boldsymbol{x}})$
$\displaystyle=\frac{1}{\bar{W}(\beta,\bar{\boldsymbol{c}};{\boldsymbol{m}})}\,\delta_{{\boldsymbol{m}}\,{\boldsymbol{m}}^{\prime}},\quad{\boldsymbol{m}},{\boldsymbol{m}}^{\prime}\in{\mathbb{N}}_{0}^{n},$
(4.9) $\displaystyle\bar{W}(\beta,\bar{\boldsymbol{c}};{\boldsymbol{m}})$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{(\beta)_{|m|}\bar{\boldsymbol{c}}^{\boldsymbol{m}}}{{\boldsymbol{m}}!}.$
(4.10)

###### Theorem 4.6

Concrete form of the polynomial in terms of ${\bf u}$ is

$\displaystyle
P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{\begin{subarray}{c}\sum_{i,j}c_{ij}\\\
(c_{ij})\in{\mathbb{M}}_{n}\end{subarray}}\frac{\prod\limits_{i=1}^{n}(-x_{i})_{\sum\limits_{j=1}^{n}c_{ij}}\prod\limits_{j=1}^{n}(-m_{j})_{\sum\limits_{i=1}^{n}c_{ij}}}{(\beta)_{\sum_{i,j}c_{ij}}}\;\frac{\prod(u_{ij})^{c_{ij}}}{\prod
c_{ij}!},$ (4.11)

in which ${\mathbb{M}}_{n}$ is the set of all $n\times n$ matrices with
entries in ${\mathbb{N}}_{0}$. This is the hypergeometric function of Aomoto-
Gelfand [1, 6] of type $(n+1,2n+2)$ and it has almost the same form as that of
the multivariate Krawtchouk polynomials reported earlier ([25].2), ([4].7),
([30].3.15), except that $-N$ is replaced by $\beta$.

It should be stressed that these results constitute the framework only. In
order to construct the multivariate Meixner polynomials explicitly, the
$n\times n$ matrix of the parameters ${\boldsymbol{u}}$ must be specified
completely. Two types of examples are presented in Iliev’s article [18]§5.

##  5 Approach via Birth and Death problem

Now I present the procedure to derive the multivariate Meixner polynomials
explicitly, i.e. to determine the parameters $\\{u_{i\,j}\\}$ explicitly,
based on the multivariate Birth and Death problem setting.

###  5.1 An overview of $n$-variate Birth and Death problem

This is a simple recapitulation of the problem setting of multivariate Birth
and Death (BD) problems reported in a previous paper [30]. For the single
variable BD problem, the well-known approach by Karlin-McGregor [21, 20] is
based on the three term recursion relations of generic orthogonal polynomials.
My method [29] employs second order difference equations which are dual to the
three term recursion relations. All of the Askey scheme polynomials of a
single discrete variable provide ‘exactly solvable BD problem’ of one
population group. For the construction of multivariate orthogonal polynomials
of Askey type, I adopt the BD problem approach after the one successful
example of the multivariate Krawtchouk polynomials [30].

The differential equation for the birth and death (BD) process of $n$ groups
of unlimited population
$\boldsymbol{x}=(x_{1},x_{2},\ldots,x_{n})\in{\mathbb{N}}_{0}$ reads

$\displaystyle\frac{\partial}{\partial t}\mathcal{P}(\boldsymbol{x};t)$
$\displaystyle=(L_{BD}\mathcal{P})(\boldsymbol{x};t)=\sum_{\boldsymbol{y}\in{\mathbb{N}}_{0}^{n}}{L_{BD}}_{\boldsymbol{x}\,\boldsymbol{y}}\mathcal{P}(\boldsymbol{y};t),\quad\mathcal{P}(\boldsymbol{x};t)\geq
0,\quad\sum_{\boldsymbol{x}\in{\mathbb{N}}_{0}^{n}}\mathcal{P}(\boldsymbol{x};t)=1,$
(5.1)
$\displaystyle=-\sum_{j=1}^{n}(B_{j}(\boldsymbol{x})+D_{j}(\boldsymbol{x}))\mathcal{P}(\boldsymbol{x};t)+\sum_{j=1}^{n}B_{j}(\boldsymbol{x}-\boldsymbol{e}_{j})\mathcal{P}(\boldsymbol{x}-\boldsymbol{e}_{j};t)$
$\displaystyle\hskip
165.02606pt+\sum_{j=1}^{n}D_{j}(\boldsymbol{x}+\boldsymbol{e}_{j})\mathcal{P}(\boldsymbol{x}+\boldsymbol{e}_{j};t),$
(5.2)

in which $\boldsymbol{e}_{j}$ is the $j$-th unit vector, $j=1,\ldots,n$. The
birth and death rates for $n$ groups $B_{j}(\boldsymbol{x})$,
$D_{j}(\boldsymbol{x})$ are all positive with the boundary conditions

$B_{j}(\boldsymbol{x})>0,\quad D_{j}(\boldsymbol{x})>0,\quad
D_{j}(\boldsymbol{x})=0\ \ \text{if}\ \ \boldsymbol{x}\in{\mathbb{N}}_{0}^{n}\
\ \text{and}\ \
\boldsymbol{x}-\boldsymbol{e}_{j}\notin{\mathbb{N}}_{0}^{n}\quad
j=1,\ldots,n.$ (5.3)

This is a typical example of the nearest neighbour interactions in $n$
dimensions. The birth and death operator $L_{BD}$, an
${\mathbb{N}}_{0}^{n}\times{\mathbb{N}}_{0}^{n}$ matrix, can be expressed
succinctly as

$\displaystyle L_{BD}$
$\displaystyle=-\sum_{j=1}^{n}\left[B_{j}(\boldsymbol{x})-B_{j}(\boldsymbol{x}-\boldsymbol{e}_{j})e^{-\partial_{j}}+D_{j}(\boldsymbol{x})-D_{j}(\boldsymbol{x}+\boldsymbol{e}_{j})e^{\partial_{j}}\right],$
(5.4)
$\displaystyle=-\sum_{j=1}^{n}\left[(1-e^{-\partial_{j}})B_{j}(\boldsymbol{x})+(1-e^{\partial_{j}})D_{j}(\boldsymbol{x})\right],$
(5.5)
$\displaystyle=-\sum_{j=1}^{n}(1-e^{-\partial_{j}})\bigl{(}B_{j}(\boldsymbol{x})-D_{j}(\boldsymbol{x}+\boldsymbol{e}_{j})\,e^{\partial_{j}}\bigr{)}.$
(5.6)

It is required that the system has a stationary distribution
$W({\boldsymbol{x}})$

$(L_{BD}W)(\boldsymbol{x})=0,\quad\sum_{\boldsymbol{x}\in{\mathbb{N}}_{0}^{n}}W(\boldsymbol{x})=1,\quad
W(\boldsymbol{x})>0,\quad\boldsymbol{x}\in{\mathbb{N}}_{0}^{n},$ (5.7)

which constrains $\\{B_{j}(\boldsymbol{x}),D_{j}(\boldsymbol{x})\\}$ severely.
The sufficient condition for the existence of the zero mode of $L_{BD}$ (5.7),
the stationary distribution $W(\boldsymbol{x})>0$, reads

$\displaystyle\bigl{(}B_{j}(\boldsymbol{x})-D_{j}(\boldsymbol{x}+\boldsymbol{e}_{j})e^{\partial_{j}}\bigr{)}W(\boldsymbol{x})=0\
\Rightarrow\frac{W(\boldsymbol{x}+\boldsymbol{e}_{j})}{W(\boldsymbol{x})}=\frac{B_{j}(\boldsymbol{x})}{D_{j}(\boldsymbol{x}+\boldsymbol{e}_{j})},\quad
j=1,\ldots,n,$ (5.8)

together with the compatibility conditions.

$\frac{B_{j}(\boldsymbol{x})}{D_{j}(\boldsymbol{x}+\boldsymbol{e}_{j})}\frac{B_{k}(\boldsymbol{x}+\boldsymbol{e}_{j})}{D_{k}(\boldsymbol{x}+\boldsymbol{e}_{j}+\boldsymbol{e}_{k})}=\frac{B_{k}(\boldsymbol{x})}{D_{k}(\boldsymbol{x}+\boldsymbol{e}_{k})}\frac{B_{j}(\boldsymbol{x}+\boldsymbol{e}_{k})}{D_{j}(\boldsymbol{x}+\boldsymbol{e}_{k}+\boldsymbol{e}_{j})},\quad
j,k=1,\ldots,n,$ (5.9)

When satisfied, these conditions determine entire $W(\boldsymbol{x})$ starting
from the origin $W(\boldsymbol{0})$.

By a similarity transformation of $L_{BD}$ in terms of
$\sqrt{W({\boldsymbol{x}})}$, a new operator (matrix) $\mathcal{H}$ is
introduced,

$\displaystyle\mathcal{H}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}-\bigl{(}\sqrt{W(\boldsymbol{x})}\bigr{)}^{-1}L_{BD}\sqrt{W(\boldsymbol{x})}.$
(5.10)
$\displaystyle=\sum_{j=1}^{n}\left[B_{j}(\boldsymbol{x})+D_{j}(\boldsymbol{x})-\sqrt{B_{j}(\boldsymbol{x})D_{j}(\boldsymbol{x}+\boldsymbol{e}_{j})}\,e^{\partial_{j}}-\sqrt{B_{j}(\boldsymbol{x}-\boldsymbol{e}_{j})D_{j}(\boldsymbol{x})}\,e^{-\partial_{j}}\right],$
(5.11) $\displaystyle\mathcal{H}_{\boldsymbol{x}\,\boldsymbol{y}}$
$\displaystyle=\sum_{j=1}^{n}\left[\bigl{(}B_{j}(\boldsymbol{x})+D_{j}(\boldsymbol{x})\bigr{)}\,\delta_{\boldsymbol{x}\,\boldsymbol{y}}-\sqrt{B_{j}(\boldsymbol{x})D_{j}(\boldsymbol{x}+\boldsymbol{e}_{j})}\,\delta_{\boldsymbol{x}+\boldsymbol{e}_{j}\,\boldsymbol{y}}\right.$
$\displaystyle\left.\hskip
147.95424pt-\sqrt{B_{j}(\boldsymbol{x}-\boldsymbol{e}_{j})D_{j}(\boldsymbol{x})}\,\delta_{\boldsymbol{x}-\boldsymbol{e}_{j}\,\boldsymbol{y}}\right].$
(5.12)

The operator $\mathcal{H}$ is a positive semi-definite real symmetric matrix
as is clear by the following factorisation,

$\displaystyle\hskip
113.81102pt\mathcal{H}=\sum_{j=1}^{n}\mathcal{A}_{j}(\boldsymbol{x})^{T}\mathcal{A}_{j}(\boldsymbol{x}),\qquad\mathcal{H}_{\boldsymbol{x}\,\boldsymbol{y}}=\mathcal{H}_{\boldsymbol{y}\,\boldsymbol{x}},$
(5.13)
$\displaystyle\mathcal{A}_{j}(\boldsymbol{x})\stackrel{{\scriptstyle\text{def}}}{{=}}\sqrt{B_{j}(\boldsymbol{x})}-e^{\partial_{j}}\sqrt{D_{j}(\boldsymbol{x})},\
\mathcal{A}_{j}(\boldsymbol{x})^{T}=\sqrt{B_{j}(\boldsymbol{x})}-\sqrt{D_{j}(\boldsymbol{x})}\,e^{-\partial_{j}},\
j=1,\ldots,n.$ (5.14)

As $W({\boldsymbol{x}})$ is the zero mode of $L_{BD}$ (5.7),
$\sqrt{W(\boldsymbol{x})}$ is the zero mode of
$\mathcal{A}_{j}(\boldsymbol{x})$ and $\mathcal{H}$

$\mathcal{A}_{j}(\boldsymbol{x})\sqrt{W(\boldsymbol{x})}=0,\quad
j=1,\ldots,n\quad\Longrightarrow\mathcal{H}\sqrt{W(\boldsymbol{x})}=0.$ (5.15)

Another operator $\widetilde{\mathcal{H}}$ is introduced by a similarity
transformation of $\mathcal{H}$ in terms of the square root of the stationary
distribution $W(\boldsymbol{x})$,

$\displaystyle\widetilde{\mathcal{H}}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\bigl{(}\sqrt{W(\boldsymbol{x})}\bigr{)}^{-1}\mathcal{H}\sqrt{W(\boldsymbol{x})}.$
(5.16)
$\displaystyle=\sum_{j=1}^{n}\left[B_{j}(\boldsymbol{x})\bigl{(}1-e^{\partial_{j}}\bigr{)}+D_{j}(\boldsymbol{x})\bigl{(}1-e^{-\partial_{j}}\bigr{)}\right],$
(5.17)

which provides the difference equations for the possible multivariate
orthogonal polynomials of Askey type. A trivial fact that a constant is the
zero mode of $\widetilde{\mathcal{H}}$ is worth mentioning

$\widetilde{\mathcal{H}}\,1=0,$ (5.18)

which corresponds to (5.15).

###  5.2 BD rates for Meixner

###### Definition 5.1

BD rates for Meixner The following Birth and Death rates are adopted for the
$n$-variate Meixner polynomials,

$B_{j}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}\beta+|x|,\qquad
D_{j}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}c_{j}^{-1}x_{j},\quad
c_{j}>0,\qquad j=1,\ldots,n,\quad\beta>0,$ (5.19)

in which the parameters $\\{c_{j}\\}$, $j=1,\ldots,n$ are asuumed to be
generic.

They trivially satisfy the compatibility conditions (5.9) for
$j,k=1,\ldots,n$,

$\displaystyle\frac{B_{j}(\boldsymbol{x})}{D_{j}(\boldsymbol{x}+\boldsymbol{e}_{j})}\frac{B_{k}(\boldsymbol{x}+\boldsymbol{e}_{j})}{D_{k}(\boldsymbol{x}+\boldsymbol{e}_{j}+\boldsymbol{e}_{k})}=\frac{c_{k}c_{k}(\beta+|x|)(\beta+|x|+1)}{(x_{j}+1)(x_{k}+1)}=\frac{B_{k}(\boldsymbol{x})}{D_{k}(\boldsymbol{x}+\boldsymbol{e}_{k})}\frac{B_{j}(\boldsymbol{x}+\boldsymbol{e}_{k})}{D_{j}(\boldsymbol{x}+\boldsymbol{e}_{k}+\boldsymbol{e}_{j})}.$

They lead to the stationary distribution introduced in Definition 4.1

$W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{(\beta)_{|x|}\boldsymbol{c}^{\boldsymbol{x}}}{\boldsymbol{x}!}(1-|c|)^{\beta},\quad\beta>0,$
(4.1)

as the sufficient conditions (5.8) are easily verified

$\displaystyle\left((\beta+|x|)-c_{j}^{-1}(x_{j}+1)e^{\partial_{j}}\right)W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})=\frac{(\beta)_{|x|+1}\boldsymbol{c}^{\boldsymbol{x}}}{\boldsymbol{x}!}(1-|c|)^{\beta}-\frac{(\beta)_{|x|+1}\boldsymbol{c}^{\boldsymbol{x}}}{\boldsymbol{x}!}(1-|c|)^{\beta}=0.$

The summability of $W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})$ due to the
summation formula

$\sum_{n=0}^{\infty}\frac{(\gamma)_{n}}{n!}z^{n}={}_{1}F_{0}\Bigl{(}\genfrac{}{}{0.0pt}{}{\gamma}{-}\Bigm{|}z\Bigr{)}=(1-|z|)^{-\gamma}$
(5.20)

limits the parameter ranges

$0<\sum_{j=1}^{n}c_{j}\equiv|c|<1.$ (5.21)

###### Definition 5.2

operator $\widetilde{\mathcal{H}}$ for $n$-variate Meixner takes a very simple
form

$\widetilde{\mathcal{H}}=(\beta+|x|)\sum_{j=1}^{n}(1-e^{\partial_{j}})+\sum_{j=1}^{n}c_{j}^{-1}x_{j}(1-e^{-\partial_{j}}).$
(5.22)

It is easy to see that the set of $n$-variate polynomials of maximal degree
$M$

$V_{M}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}\text{Span}\\{{\boldsymbol{x}}^{\boldsymbol{m}}|0\leq|m|\leq
M\\},\quad|m|\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{j=1}^{n}m_{j},$
(5.23)

is invariant under $\widetilde{\mathcal{H}}$

$\widetilde{\mathcal{H}}V_{M}({\boldsymbol{x}})\subseteq
V_{M}({\boldsymbol{x}}),$ (5.24)

and $\widetilde{\mathcal{H}}$ has eigenpolynomials in each
$V_{M}({\boldsymbol{x}})$. For generic values of the parameters $\\{c_{j}\\}$,
these polynomials are orthogonal with each other due to the real symmetry of
the operator $\mathcal{H}$ (5.13).

###  5.3 Meixner polynomials are orthogonal with each other

It is easy to determine $n$ degree 1 eigenpolynomials of
$\widetilde{\mathcal{H}}$ (5.22) with unknown coefficients $\\{a_{i}\\}$ and
unit constant part,

$\displaystyle
P_{|m|=1}(\boldsymbol{x})=1+\sum_{i=1}^{n}a_{i}x_{i},\quad\widetilde{\mathcal{H}}P_{|m|=1}(\boldsymbol{x})=\lambda
P_{|m|=1}(\boldsymbol{x}),$ (5.25)
$\displaystyle\Rightarrow-(\beta+|x|)\sum_{i=1}^{n}a_{i}+\sum_{i=1}^{n}c_{i}^{-1}a_{i}x_{i}=\lambda\Bigl{(}1+\sum_{i=1}^{n}a_{i}x_{i}\Bigr{)}.$

By equating the coefficients of $x_{i}$ and 1, an eigenvalue equations of
$\\{a_{i}\\}$ are obtained,

$\displaystyle-\sum_{k=1}^{n}a_{k}+c_{i}^{-1}a_{i}$ $\displaystyle=\lambda
a_{i},\qquad-\beta\sum_{k=1}^{n}a_{k}=\lambda,$ (5.26)
$\displaystyle\Longrightarrow\ a_{i}=\frac{1}{\beta}\frac{\lambda}{\lambda-
c_{i}^{-1}},\quad i=1,\ldots,n,.$ (5.27)

Here $\lambda$ is the root of a degree $n$ characteristic polynomial
$\mathcal{F}(\lambda)$ of an $n\times n$ matrix $F({\boldsymbol{c}})$
depending on $\\{c_{i}\\}$,

$0=\mathcal{F}(\lambda)\stackrel{{\scriptstyle\text{def}}}{{=}}\text{Det}\bigl{(}\lambda
I_{n}-F({\boldsymbol{c}})\bigr{)},\quad
F({\boldsymbol{c}})_{i\,j}\stackrel{{\scriptstyle\text{def}}}{{=}}-1+c_{i}^{-1}\delta_{i\,j}.$
(5.28)

For each eigenvalue $\lambda_{j}$, which is positive by construction, the
unknown coefficients $\\{a_{i}\\}$ are determined,

$a_{i,j}=\frac{\lambda_{j}}{\beta(\lambda_{j}-c_{i}^{-1})},\quad
i,j=1,\ldots,n,$

and it satisfies the relation

$\sum_{i=1}^{n}\frac{1}{\lambda_{j}-c_{i}^{-1}}\equiv\sum_{i=1}^{n}\frac{c_{i}}{c_{i}\lambda_{j}-1}=-1,\quad
j=1,\ldots,n.$ (5.29)

Let us tentatively identify the above $j$-th solution as
$\boldsymbol{m}=\boldsymbol{e}_{j}$ solution

$P_{\boldsymbol{e}_{j}}(\boldsymbol{x})=1+\frac{1}{\beta}\sum_{i=1}^{n}\frac{\lambda_{j}}{\lambda_{j}-c_{i}^{-1}}x_{i},\quad
j=1,\ldots,n.$ (5.30)

By comparing these polynomials with the general hypergeometric functions [1,
6] in Theorem 4.6

$\displaystyle
P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{\begin{subarray}{c}\sum_{i,j}c_{ij}\\\
(c_{ij})\in{\mathbb{M}}_{n}\end{subarray}}\frac{\prod\limits_{i=1}^{n}(-x_{i})_{\sum\limits_{j=1}^{n}c_{ij}}\prod\limits_{j=1}^{n}(-m_{j})_{\sum\limits_{i=1}^{n}c_{ij}}}{(\beta)_{\sum_{i,j}c_{ij}}}\;\frac{\prod(u_{ij})^{c_{ij}}}{\prod
c_{ij}!},$ (4.11)

the system parameters $\\{u_{i\,j}\\}$ are completely identified

$P_{\boldsymbol{e}_{j}}(\beta,{\boldsymbol{u}};\boldsymbol{x})=1+\frac{1}{\beta}\sum_{i=1}^{n}u_{i\,j}x_{i},\qquad
u_{i\,j}=\frac{\lambda_{j}}{\lambda_{j}-c_{i}^{-1}},\quad i,j=1,\ldots,n.$
(5.31)

The $n+1$ eigenvectors $\sqrt{W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})}$,
$\\{\sqrt{W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})}P_{\boldsymbol{e}_{j}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\\}$,
$j=1,\ldots,n$ of the real symmetric matrix $\mathcal{H}$ (5.11) are
orthogonal with each other for generic parameters $\\{c_{i}\\}$,

$\displaystyle\sum_{{\boldsymbol{x}}\in{\mathbb{N}}_{0}^{n}}W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})P_{\boldsymbol{e}_{j}}(\beta,{\boldsymbol{u}};\boldsymbol{x})$
$\displaystyle=0,\qquad\qquad\quad j=1,\ldots,n,$ (5.32)
$\displaystyle\sum_{{\boldsymbol{x}}\in{\mathbb{N}}_{0}^{n}}W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})P_{\boldsymbol{e}_{j}}(\beta,{\boldsymbol{u}};\boldsymbol{x})P_{\boldsymbol{e}_{k}}(\beta,{\boldsymbol{u}};\boldsymbol{x})$
$\displaystyle=0,\quad j\neq k,\quad j,k=1,\ldots,n.$ (5.33)

By using the summation formulas

$\displaystyle\sum_{{\boldsymbol{x}}\in{\mathbb{N}}_{0}^{n}}W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})x_{j}$
$\displaystyle=\frac{c_{j}\beta}{1-|c|},\qquad\qquad\qquad\qquad\quad\
j=1,\ldots,n,$
$\displaystyle\sum_{{\boldsymbol{x}}\in{\mathbb{N}}_{0}^{n}}W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})x_{j}x_{k}$
$\displaystyle=\frac{\beta(\beta+1)c_{j}c_{k}}{(1-|c|)^{2}}+\frac{\beta
c_{j}}{1-|c|}\,\delta_{j\,k},\quad j,k=1,\ldots,n,$

they are reduced to the same expressions as the constraining conditions (4.3)
and (4.4) of the $n^{2}$ real parameters $\boldsymbol{u}=\\{u_{i\,j}\\}$,
$i,j=1,\ldots,n$ in Definition 4.1,

$\displaystyle\sum_{i=1}^{n}c_{i}u_{i\,j}=|c|-1,\quad j=1,\ldots,n,$ (4.3)
$\displaystyle\sum_{i=1}^{n}c_{i}u_{i\,j}u_{i\,k}=|c|-1,\quad j\neq k,\quad
j,k=1,\ldots,n.$ (4.4)

The norm of $P_{\boldsymbol{e}_{j}}(\beta,{\boldsymbol{u}};\boldsymbol{x})$
can be calculated similarly

$\sum_{{\boldsymbol{x}}\in{\mathbb{N}}_{0}^{n}}W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})P_{\boldsymbol{e}_{j}}(\beta,{\boldsymbol{u}};\boldsymbol{x})^{2}=\frac{1}{\beta}\frac{1-|c|+\sum_{i=1}^{n}c_{i}u_{i\,j}^{2}}{1-|c|}=\frac{1}{\beta\bar{c}_{j}},\qquad
j=1,\ldots,n.$ (4.6)

These lead to the following

###### Theorem 5.3

Orthogonality of $P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})$
The explicitly assembled parameters
$\\{u_{i\,j}=\lambda_{j}/(\lambda_{j}-c_{i}^{-1})\\}$, $i,j=1,\ldots,n$ (5.31)
satisfy all the constraints required for the definition of the generating
function $G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})$ (4.7) in
Definition 4.3. Definition 4.4 and Theorem 4.5, 4.6 assure that the
polynomials $P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})$ (4.11)
with the above ${\boldsymbol{u}}=\\{u_{i\,j}\\}$ satisfy the orthogonality
relation

$\displaystyle\sum_{{\boldsymbol{x}}\in{\mathbb{N}}_{0}^{n}}W(\beta,{\boldsymbol{c}};{\boldsymbol{x}})P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};{\boldsymbol{x}})P_{{\boldsymbol{m}}^{\prime}}(\beta,{\boldsymbol{u}};{\boldsymbol{x}})$
$\displaystyle=\frac{1}{\bar{W}(\beta,\bar{\boldsymbol{c}};{\boldsymbol{m}})}\,\delta_{{\boldsymbol{m}}\,{\boldsymbol{m}}^{\prime}},\quad{\boldsymbol{m}},{\boldsymbol{m}}^{\prime}\in{\mathbb{N}}_{0}^{n},$
(4.9) $\displaystyle\bar{W}(\beta,\bar{\boldsymbol{c}};{\boldsymbol{m}})$
$\displaystyle=\frac{(\beta)_{|m|}\bar{\boldsymbol{c}}^{\boldsymbol{m}}}{{\boldsymbol{m}}!},$
(4.10)

with the $\\{\bar{c}_{j}>0\\}$, $j=1,\ldots,n$ determined explicitly by (4.6).

###  5.4 All
$\\{P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};{\boldsymbol{x}})\\}$ are
eigenvectors of $\widetilde{\mathcal{H}}$

Verifying that all higher degree ones
$\\{P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\\}$ (4.11) are
also the eigenpolynomials of $\widetilde{\mathcal{H}}$ (5.22) is the next
task. The explicit forms of the eigenvalues $\mathcal{E}(\boldsymbol{m})$ are
necessary. Since $P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};{\boldsymbol{x}})$
has the form

$P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})=1+\frac{1}{\beta}\sum_{i,j=1}^{n}x_{i}m_{j}u_{i\,j}+\text{higher
degrees},$ (5.34)

$\widetilde{\mathcal{H}}$ acting on the higher degrees produces only the terms
of linear and higher degrees. The only constant part of
$\widetilde{\mathcal{H}}P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})$
comes from $\beta\sum_{i=1}^{n}(1-e^{\partial_{j}})$ acting on the linear
part,

$\beta\sum_{k=1}^{n}(1-e^{\partial_{k}})\left\\{\frac{1}{\beta}\sum_{i,j=1}^{n}x_{i}m_{j}u_{i\,j}\right\\}=-\sum_{i\,j}m_{j}u_{i\,j}=-\sum_{j=1}^{n}m_{j}\lambda_{j}\sum_{i=1}^{n}\frac{1}{\lambda_{j}-c_{i}^{-1}}=\sum_{j=1}^{n}m_{j}\lambda_{j},$

in which (5.29) is used. After applying $(1-e^{\partial_{i}})$, all higher
degree terms vanish at the origin ${\boldsymbol{x}}={\boldsymbol{0}}$ as they
consist of terms like $(x_{i})_{k}(x_{j})_{l}$, $k+l\geq 2$, the typical
structure of the hypergeometric functions. This leads to the following

###### Proposition 5.4

If $P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})$ (4.11) is an
eigenpolynomial of $\widetilde{\mathcal{H}}$ (5.22), it has a linear spectrum

$\mathcal{E}(\boldsymbol{m})\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{j=1}^{n}m_{j}\lambda_{j}.$

The remaining task is to demonstrate

$\widetilde{\mathcal{H}}P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})=\Bigl{(}\sum_{k=1}^{n}m_{k}\lambda_{k}\Bigr{)}P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x}),\qquad\boldsymbol{x},\boldsymbol{m}\in{\mathbb{N}}_{0}^{n}.$
(5.35)

After the example of the single variable Meixner polynomials in §2, the
following $n$-variables generalisation of the generating function formula
(2.7)

$\displaystyle\widetilde{\mathcal{H}}G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})$
$\displaystyle=\left(\sum_{k=1}^{n}\lambda_{k}t_{k}\frac{\partial}{\partial
t_{k}}\right)G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t}),$ (5.36)
$\displaystyle
G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})\stackrel{{\scriptstyle\text{def}}}{{=}}(1-|t|)^{-\beta-|x|}\prod_{i=1}^{n}\left(1-\sum_{j=1}^{n}b_{i\,j}t_{j}\right)^{x_{i}},$
(4.7)

will lead to (5.35) above.

The action of the two parts of $\widetilde{\mathcal{H}}$

$(\beta+|x|)\sum_{i=1}^{n}(1-e^{\partial_{i}}),\qquad\sum_{i=1}^{n}c_{i}^{-1}x_{i}(1-e^{-\partial_{i}}),$

on $G(u,\boldsymbol{x},t)$ is evaluated separately. The first part gives

$\displaystyle(\beta+|x|)\sum_{i=1}^{n}(1-e^{\partial_{i}})G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})$
$\displaystyle\qquad=(\beta+|x|)\sum_{i=1}^{n}\left(G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})-\frac{1-\sum_{j=1}^{n}b_{i\,j}t_{j}}{1-|t|}G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})\right)$
$\displaystyle\qquad=(\beta+|x|)\frac{1}{1-|t|}\sum_{i,j=1}^{n}(b_{i\,j}-1)t_{j}G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})$
$\displaystyle\qquad=(\beta+|x|)\sum_{j=1}^{n}\frac{\lambda_{j}t_{j}}{1-|t|}G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t}),\hskip
128.0374pt(*)$

in which (4.5), (5.31) and (5.29) are used to obtain

$\sum_{i=1}^{n}\sum_{j=1}^{n}(b_{i\,j}-1)t_{j}=-\sum_{i=1}^{n}\sum_{j=1}^{n}u_{i\,j}t_{j}=-\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{\lambda_{j}}{\lambda_{j}-c_{i}^{-1}}t_{j}=\sum_{j=1}^{n}\lambda_{j}t_{j}.$

The second part gives

$\displaystyle\sum_{i=1}^{n}c_{i}^{-1}x_{i}(1-e^{\partial_{i}})G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})$
$\displaystyle\qquad=\sum_{i=1}^{n}c_{i}^{-1}x_{i}\left(G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})-\frac{1-|t|}{1-\sum_{j=1}^{n}b_{i\,j}t_{j}}G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})\right)$
$\displaystyle\qquad=\sum_{i=1}^{n}c_{i}^{-1}x_{i}\frac{|t|-\sum_{j=1}^{n}b_{i\,j}t_{j}}{1-\sum_{j=1}^{n}b_{i\,j}t_{j}}G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})$
$\displaystyle\qquad=\sum_{i,j=1}^{n}\frac{c_{i}^{-1}x_{i}u_{i\,j}t_{j}}{1-\sum_{j=1}^{n}b_{i\,j}t_{j}}G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t}).\hskip
128.0374pt(**)$

Now r.h.s. of (5.36) reads

$\displaystyle(\beta+|x|)\sum_{k=1}^{n}\frac{\lambda_{k}t_{k}}{1-|t|}G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t})-\sum_{i,k=1}^{n}\frac{\lambda_{k}t_{k}x_{i}b_{i\,k}}{1-\sum_{j=1}^{n}b_{i\,j}t_{j}}G(\beta,\boldsymbol{u},\boldsymbol{x};\boldsymbol{t}).$

Here the following equality holds

$\lambda_{k}x_{i}b_{i\,k}=\lambda_{k}x_{i}\frac{-c_{i}^{-1}}{\lambda_{k}-c_{i}^{-1}}=-c_{i}^{-1}x_{i}\frac{\lambda_{k}}{\lambda_{k}-c_{i}^{-1}}=c_{i}^{-1}x_{i}u_{i\,k}.$

The r.h.s. are qual to $(*)+(**)$ and this concludes the proof, which leads to
the following

###### Theorem 5.5

The multivariate Meixner polynomials
$\\{P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\\}$ (4.11)
constitute the complete set of eigenpolynomils of the difference operator
$\widetilde{\mathcal{H}}$ (5.22) which is derived from the $n$-variate Birth
and Death process (5.19). The orthogonality (4.9) is the consequence of the
self-adjointness of $\mathcal{H}$ (5.12) and the generic parameters.

###### Proposition 5.6

$\mathfrak{S}_{n}$ Symmetry The Meixner polynomials
$\\{P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\\}$ (4.11) is
invariant under the symmetric group $\mathfrak{S}_{n}$, due to the
arbitrariness of the ordering of $n$ roots $\\{\lambda_{j}\\}$ of the
characteristic equation (5.28) in the parameters $u_{i\,j}$ (5.31).

###  5.5 Solution of the BD problem

In terms of the norm formula (4.9), the following set of orthonormal vectors
on the Hilbert space ${\mathbb{N}}_{0}^{n}$ are defined

$\displaystyle\quad\sum_{\boldsymbol{x}\in{\mathbb{N}}_{0}^{n}}\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{m}^{\prime}}(\beta,{\boldsymbol{u}};\boldsymbol{x})=\delta_{\boldsymbol{m}\,\boldsymbol{m}^{\prime}},\qquad\quad\boldsymbol{m},\boldsymbol{m}^{\prime}\in{\mathbb{N}}_{0}^{n},,$
(5.37)
$\displaystyle\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\sqrt{W(\beta,{\boldsymbol{c}};\boldsymbol{x})}P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\sqrt{\bar{W}(\beta,\bar{\boldsymbol{c}};\boldsymbol{m})},\qquad\boldsymbol{x},\boldsymbol{m}\in{\mathbb{N}}_{0}^{n},$
(5.38)
$\displaystyle\bar{W}(\beta,\bar{\boldsymbol{c}};{\boldsymbol{m}})=\frac{(\beta)_{|m|}\bar{\boldsymbol{c}}^{\boldsymbol{m}}}{{\boldsymbol{m}}!},\quad\sum_{\boldsymbol{m}\in{\mathbb{N}}_{0}^{n}}\bar{W}(\beta,\bar{\boldsymbol{c}};{\boldsymbol{m}})=(1-|\bar{c}|)^{-\beta}.$
(5.39)

They define an orthogonal matrix $\mathcal{S}$ on ${\mathbb{N}}_{0}^{n}$,

$\displaystyle\hskip
113.81102pt\mathcal{S}_{\boldsymbol{x}\,\boldsymbol{m}}\stackrel{{\scriptstyle\text{def}}}{{=}}\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x}),$
$\displaystyle\bigl{(}\mathcal{S}^{T}\mathcal{S}\bigr{)}_{\boldsymbol{m}\,\boldsymbol{m}^{\prime}}=\sum_{\boldsymbol{x}\in{\mathbb{N}}_{0}^{n}}\bigl{(}\mathcal{S}^{T}\bigr{)}_{\boldsymbol{m}\,\boldsymbol{x}}\mathcal{S}_{\boldsymbol{x}\,\boldsymbol{m}^{\prime}}=\sum_{\boldsymbol{x}\in{\mathbb{N}}_{0}^{n}}\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{m}^{\prime}}(\beta,{\boldsymbol{u}};\boldsymbol{x})=\delta_{\boldsymbol{m}\,\boldsymbol{m}^{\prime}},$
$\displaystyle\Rightarrow\delta_{\boldsymbol{x}\,\boldsymbol{y}}=\bigl{(}\mathcal{S}\mathcal{S}^{T})_{\boldsymbol{x}\,\boldsymbol{y}}=\sum_{\boldsymbol{m}\in{\mathbb{N}}_{0}^{n}}\mathcal{S}_{\boldsymbol{x}\,\boldsymbol{m}}\bigl{(}\mathcal{S}^{T}\bigr{)}_{\boldsymbol{m}\,\boldsymbol{y}}=\sum_{\boldsymbol{m}\in{\mathbb{N}}_{0}^{n}}\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{y}).$
(5.40)

This means that
$\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})$ defines
dual polynomials in $\boldsymbol{m}$ indexed by $\boldsymbol{x}$. The
discussion of dual polynomials is essentially the same as that for the
multivariate Krawtchouk polynomials in [30] and it will not be repeated here.

In terms of the orthonormal eigenvectors
$\\{\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\\}$ of
$\mathcal{H}$ (5.11) and $L_{BD}$ (5.4) the transition probability
$\mathcal{T}({\boldsymbol{x}},{\boldsymbol{y}};t)$ of the Birth and Death
equation (5.1) is expressed neatly. The transition probability
$\mathcal{T}({\boldsymbol{x}},{\boldsymbol{y}};t)$ is the solution of the BD
equation (5.1) with the initial condition

$\mathcal{P}(\boldsymbol{x};0)=\delta_{\boldsymbol{x}\,\boldsymbol{y}},$
$\displaystyle\mathcal{T}({\boldsymbol{x}},{\boldsymbol{y}};t)$
$\displaystyle=\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{y})^{-1}\sum_{{\boldsymbol{m}}\in{\mathbb{N}}_{0}^{n}}\,e^{-\mathcal{E}(\boldsymbol{m})t}\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{y}),\quad
t>0.$ (5.41)

It is trivial to verify the initial condition by (5.40). By using

$\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{x})=\sqrt{W(\beta,{\boldsymbol{c}};\boldsymbol{x})}\sqrt{\bar{W}(\beta,\bar{\boldsymbol{c}};\boldsymbol{m})},$

(5.41) is reduced to

$\displaystyle\mathcal{T}({\boldsymbol{x}},{\boldsymbol{y}};t)$
$\displaystyle=W(\beta,{\boldsymbol{c}};\boldsymbol{x})\bar{W}(\beta,\bar{\boldsymbol{c}};\boldsymbol{m})\sum_{{\boldsymbol{m}}\in{\mathbb{N}}_{0}^{n}}\,e^{-\mathcal{E}(\boldsymbol{m})t}P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{y}),\quad
t>0,$ $\displaystyle\frac{\partial}{\partial
t}\mathcal{T}({\boldsymbol{x}},{\boldsymbol{y}};t)$
$\displaystyle=-W(\beta,{\boldsymbol{c}};\boldsymbol{x})\bar{W}(\beta,\bar{\boldsymbol{c}};\boldsymbol{m})\sum_{{\boldsymbol{m}}\in{\mathbb{N}}_{0}^{n}}\,\mathcal{E}(\boldsymbol{m})\,e^{-\mathcal{E}(\boldsymbol{m})t}P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{y}),\quad
t>0.$

Since
$P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})W(\beta,{\boldsymbol{c}};\boldsymbol{x})$
is the eigenvector of $L_{BD}$ with the eigenvalue
$-\mathcal{E}(\boldsymbol{m})$,

$L_{BD}P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})W(\beta,{\boldsymbol{c}};\boldsymbol{x})=-\mathcal{E}(\boldsymbol{m})P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})W(\beta,{\boldsymbol{c}};\boldsymbol{x}),$

the BD equation is satisfied

$\frac{\partial}{\partial
t}\mathcal{T}({\boldsymbol{x}},{\boldsymbol{y}};t)=L_{BD}\mathcal{T}({\boldsymbol{x}},{\boldsymbol{y}};t),\quad
t>0.$

The transition probability (5.41) ptrovides a simple example of Chapman-
Kolmogorov equation in ${\mathbf{N}}_{0}^{n}$,

$\mathcal{T}({\boldsymbol{x}},{\boldsymbol{y}};t+t^{\prime})=\sum_{{\boldsymbol{z}}\in{\mathbb{N}}_{0}^{n}}\mathcal{T}({\boldsymbol{x}},{\boldsymbol{z}};t)\mathcal{T}({\boldsymbol{z}},{\boldsymbol{y}};t^{\prime}).$
(5.42)

The r.h.s. is

$\displaystyle\sum_{{\boldsymbol{z}}\in{\mathbb{N}}_{0}^{n}}\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{z})^{-1}\sum_{{\boldsymbol{m}}\in{\mathbb{N}}_{0}^{n}}\,e^{-\mathcal{E}(\boldsymbol{m})t}\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{z})$
$\displaystyle\qquad\times\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{z})\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{y})^{-1}\sum_{{\boldsymbol{m}}^{\prime}\in{\mathbb{N}}_{0}^{n}}\,e^{-\mathcal{E}(\boldsymbol{m}^{\prime})t^{\prime}}\hat{\phi}_{\boldsymbol{m}^{\prime}}(\beta,{\boldsymbol{u}};\boldsymbol{z})\hat{\phi}_{\boldsymbol{m}^{\prime}}(\beta,{\boldsymbol{u}};\boldsymbol{y})$
$\displaystyle=\sum_{{\boldsymbol{z}}\in{\mathbb{N}}_{0}^{n}}\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{y})^{-1}\sum_{{\boldsymbol{m}}\in{\mathbb{N}}_{0}^{n}}\,e^{-\mathcal{E}(\boldsymbol{m})t}\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{z})$
$\displaystyle\qquad\times\sum_{{\boldsymbol{m}}^{\prime}\in{\mathbb{N}}_{0}^{n}}\,e^{-\mathcal{E}(\boldsymbol{m}^{\prime})t^{\prime}}\hat{\phi}_{\boldsymbol{m}^{\prime}}(\beta,{\boldsymbol{u}};\boldsymbol{z})\hat{\phi}_{\boldsymbol{m}^{\prime}}(\beta,{\boldsymbol{u}};\boldsymbol{y}).$

The summation over ${\boldsymbol{z}}$,
$\sum_{{\boldsymbol{z}}\in{\mathbb{N}}_{0}^{n}}\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{z})\hat{\phi}_{\boldsymbol{m}^{\prime}}(\beta,{\boldsymbol{u}};\boldsymbol{z})=\delta_{{\boldsymbol{m}}\,{\boldsymbol{m}}^{\prime}}$,
gives

$\displaystyle r.h.s.$
$\displaystyle=\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{0}}(\beta,{\boldsymbol{u}};\boldsymbol{y})^{-1}\sum_{{\boldsymbol{m}}\in{\mathbb{N}}_{0}^{n}}\,e^{-\mathcal{E}(\boldsymbol{m})(t+t^{\prime})}\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\hat{\phi}_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{y})$
$\displaystyle=l.h.s..$

###  5.6 Exceptional Cases

So far the parameter values are assumed to be generic. But obviously at
certain parameter settings, the above hypergeometric formula (4.11) for the
polynomials $\\{P_{\boldsymbol{m}}(\beta,{\boldsymbol{u}};\boldsymbol{x})\\}$
could go wrong. By construction, the eigenvalues $\\{\lambda_{i}\\}$ are
positive and $\\{c_{i}\\}$ are also positive. Therefore, if the situation
$\lambda_{j}=c_{i}^{-1}$ happens at some parameter setting, it leads to the
breakdown of the generic theory as
$u_{i\,j}=\frac{\lambda_{j}}{\lambda_{j}-c_{i}^{-1}}$ (5.31) is ill-defined.

####  5.6.1 $n=2$ Case

The situation is most clearly seen when $n=2$. In this case the two
eigenvalues are the roots of the quadratic equation

$\displaystyle\lambda^{2}+(2-c_{1}^{-1}-c_{2}^{-1})\lambda+(1-c_{1}^{-1})(1-c_{2}^{-1})-1=0,$
$\displaystyle\lambda_{1}$
$\displaystyle=\frac{1}{2}(-2+c_{1}^{-1}+c_{2}^{-1}-\Delta),\quad\lambda_{2}=\frac{1}{2}(-2+c_{1}^{-1}+c_{2}^{-1}+\Delta),$
$\displaystyle\Delta^{2}=4+(1/c_{1}-1/c_{2})^{2}.$

When $c_{1}=c_{2}=c$, the eigenvalues are rational,

$c_{1}=c_{2}=c\ \Longrightarrow\lambda_{1}=-2+1/c>0,\quad\lambda_{2}=1/c,$

and the singular situation occurs. That is, the general formula (4.11) fails.
In this case, the degree 1 solution of $\widetilde{\mathcal{H}}$ (5.22)
corresponding to the eigenvalue $\lambda_{2}=1/c$ is

$P_{\boldsymbol{e}_{1}}(\boldsymbol{x})=const\times\bigl{(}x_{1}-x_{2}).$

That is the constant part is vanishing and the assumption that degree one
solutions have unit constant part (5.25) simply fails. The situation is
similar for general $n$ as stated by the following

###### Theorem 5.7

When some of the parameters $\\{c_{i}\\}$ coincide, the hypergeometric formula
(4.11) for the $n$-variate Meixner polynomials does not apply. But the
solutions of the difference equations $\widetilde{\mathcal{H}}$ (5.22) still
constitute the $n$-variate orthogonal polynomials.

This is rather easy to see. If $c_{j}=c_{k}=c$, the matrix
$c^{-1}I_{n}-F({\boldsymbol{c}})$ (5.28) has the $j$-th column
$-(1,1,\ldots,1)^{T}$ and $k$-th column $-(1,1,\ldots,1)^{T}$, thus the
characteristic polynomial $\mathcal{F}(\lambda)$ vanishes at $\lambda=c^{-1}$.
In this case it is easy to show

$\widetilde{\mathcal{H}}(x_{j}-x_{k})=\frac{1}{c}(x_{j}-x_{k}).$

When there exist $k$ identical $c_{i}$’s, $\mathcal{F}(\lambda)$ has a factor
$(\lambda-c_{i}^{-1})^{k-1}$.

###### Theorem 5.8

Distinct parameters $\\{c_{j}\\}$ are necessary
All the parameters $\\{c_{j}\\}$ must be distinct for the hypergeometric
formula (4.11) for the $n$-variate Meixner polynomials to hold.

###### Remark 5.9

It is a big challenge to derive a general formula of $n$-variate Meixner
polynomials including all these exceptional cases.

## References

  * [1] Aomoto K., Kita M., Hypergeometric funtions, Springer, Berlin, 1994 (in Japanese).
  * [2] Andrews G. E., Askey R. and Roy R., Special Functions, Encyclopedia of mathematics and its applications, Cambridge Univ. Press, Cambridge, (1999).
  * [3] Cooper R. D., Hoare M. R. and Rahman M., “Stochastic Processes and Special Functions: On the Probabilistic Origin of Some Positive Kernels Associated with Classical Orthogonal Polynomials,” J. Math. Anal. Appl. 61 (1977) 262-291.
  * [4] Diaconis P. and Griffiths R.C., “An introduction to multivariate Krawtchouk polynomials and their applications,” J. Stat. Planning and Inference 154 (2014) 39-53, arXiv:1309.0112[math.PR].
  * [5] Feller W., An Introduction to Probability Theory and its Applications, I, (2nd ed.), Wiley, New York, (1966).
  * [6] Gelfand I. M., “General theory of hypergeometric functions,” Sov. Math. Dokl. 33 (1986) 573-577.
  * [7] Genest V. X., Miki H., Vinet L. and Zhedanov A., “Spin lattices, state transfer and bivariate Krawtchouk polynomials,” J. Phys. A: Math. Theor. 46 (2013) 505203, arXiv:1410.4703[math-ph].
  * [8] Griffiths R. C., “Orthogonal polynomials on the multinomial distribution,” Austral. J. Statist. 13 (1971) 27–35.
  * [9] Griffiths R. “Lancaster distributions and Markov chains with multivariate Poisson-Charlier, Meixner and Hermite-Chebycheff polynomial eigenfunctions,” J. Approxi. Theo. 207 (2016) 139-164, arXiv:1412.3931[math.PR].
  * [10] Griffiths R. C., “Multivariate Krawtchouk Polynomials and Composition Birth and Death Processes,” Symmetry 8 (2016) 33 19pp. arXiv:1603.00196[math.PR].
  * [11] Grünbaum F. A., “The Rahman polynomials are bispectral,” SIGMA 3 (2007) 065, 11pp, arXiv:0705.0468[math.CA].
  * [12] Grünbaum F. A. and Rahman M., “On a family of 2-variable orthogonal Krawtchouk polynomials,” SIGMA 6 (2010) 090, 12 pages, arXiv:1007.4327[math.CA].
  * [13] Grünbaum F. A. and Rahman M., “A System of multivariable Krawtchouk polynomials and a probabilistic application,” SIGMA 7 (2011) 119, 17pp, arXiv:1106.1835[math.PR].
  * [14] Hoare M. R. and Rahman M., “A probabilistic origin for a new class of bivariate polynomials,” SIGMA 4 (2008) 089, 18 pages, arXiv:0812.3879[math.CA].
  * [15] Iliev P. and Xu Y., “Discrete orthogonal polynomials and difference equations of several variables,” Adv. Math. 212 (2007) 1-36, arXiv:math.CA/0508039.
  * [16] Iliev P. and Terwilliger P., “The Rahman polynomials and the Lie algebra $sl_{3}(C)$,” Trans. Amer. Math. Soc. 364 (2012) 4225–4238, arXiv:1006.5062[math.RT].
  * [17] Iliev P., “A Lie theoretic interpretation of multivariate hypergeometric polynomials,” Compositio Math. 148 (2012) 991-1002, arXiv:1101.1683[math.RT].
  * [18] Iliev P., “Meixner polynomials in several variables satisfying bispectral difference equations,” Adv. in Appl. Math. 49 (2012) 15-23, arXiv:1112.5589[math.CA]
  * [19] Iliev P., “Gaudin model for the multinomial distribution,” arXiv:2303.08206[math-ph].
  * [20] Ismail M. E. H., Classical and quantum orthogonal polynomials in one variable, Encyclopedia of mathematics and its applications, Cambridge, (2005).
  * [21] Karlin S. and McGregor J. L., “The differential equations of birth-and-death processes,” Trans. Amer. Math. Soc. 85 (1957) 489-546; “Linear growth, birth-and-death processes,” J. Math. Mech. 7 (1958) 643–662; “Ehrenfest urn models,” J. Appl. Prob. 19 (1965) 477–487.
  * [22] Khare K. and Zhou H., “Rates of convergence of some multivariate Markov chains with polynomial eigenfucntions,” Ann. Appl. Probab. 19 (2009) 737-777, arXiv:0906.4242[math.PR].
  * [23] Koekoek R., Lesky P. A. and Swarttouw R. F., Hypergeometric orthogonal polynomials and their $q$-analogues, Springer Monographs in Mathematics, Springer-Verlag, Berlin, (2010).
  * [24] Mizukawa H., “Zonal spherical functions on the complex reflection groups and $(m+1,n+1)$-hypergeometric functions,” Adv. Math. 184 (2004) 1–17.
  * [25] Mizukawa H., “Orthogonal relations for multivariate Krawtchouk polynomials,” SIGMA 7 (2011) 017 5pp, arXiv:1009.1203[math.CO].
  * [26] Mizukawa H. and Tanaka H., $(n+1,m+1)$-hypergeometric functions associated to character algebras, Proc. Amer. Math. Soc. 132 (2004) 2613–2618.
  * [27] Odake S. and Sasaki R., “Orthogonal Polynomials from Hermitian Matrices,” J. Math. Phys. 49 (2008) 053503 (43 pp), arXiv:0712.4106[math.CA].
  * [28] Odake S. and Sasaki R., “Orthogonal Polynomials from Hermitian Matrices II,” J. Math. Phys. 59 (2018) 013504 (42pp) arXiv:1604.00714[math.CA].
  * [29] Sasaki R., “Exactly Solvable Birth and Death Processes,” J. Math. Phys. 50 (2009) 103509 (18 pp), arXiv:0903.3097[math-ph].
  * [30] Sasaki R., “Multivariate Krawtchouk polynomials as Birth and Death polynomials,” arXiv:2305.08581v2[math.CA]
  * [31] Tratnik M. V., “Some multivariable orthogonal polynomials of the Askey tableau-discrete families,” J. Math. Phys. 32 (1991), 2337-2342.
  * [32] Xu Y., “Hahn, Jacobi, and Krawtchouk polynomials of several variables,” Journal of Approximation Theory 195 (2015) 19-42, arXiv:1309.1510[math.CA].
  * [33] Zhedanov A., “9j-symbols of the oscillator algebra and Krawtchouk polynomials in two variables,” J. Phys. A: Math. Gen. 30 (1997) 8337-8353.