Multidimensional Toda Lattices: Continuous and Discrete Time

Multidimensional Toda Lattices:
Continuous and Discrete Time\star\star$\star$This paper is a contribution to
the Special Issue on Orthogonal Polynomials, Special Functions and
Applications. The full collection is available at
http://www.emis.de/journals/SIGMA/OPSFA2015.html

Alexander I. APTEKAREV ${}^{{\dagger}^{1}}$, Maxim DEREVYAGIN
${}^{{\dagger}^{2}}$, Hiroshi MIKI ${}^{{\dagger}^{3}}$
and Walter VAN ASSCHE ${}^{{\dagger}^{4}}$

A.I. Aptekarev, M. Derevyagin, H. Miki and W. Van Assche

${}^{{\dagger}^{1}}$ Keldysh Institute for Applied Mathematics, Russian
Academy of Sciences,
${}^{{\dagger}^{1}}$ Miusskaya pl. 4, 125047 Moscow, Russia<EMAIL_ADDRESS>
${}^{{\dagger}^{2}}$ University of Mississippi, Department of Mathematics,
${}^{{\dagger}^{2}}$ Hume Hall 305, P. O. Box 1848, University, MS 38677-1848,
USA<EMAIL_ADDRESS>
${}^{{\dagger}^{3}}$ Doshisha University, Department of Electronics, Faculty
of Science and Engineering,
${}^{{\dagger}^{3}}$ Kyotanabe city, Kyoto 610 0394, Japan
<EMAIL_ADDRESS>
${}^{{\dagger}^{4}}$ KU Leuven, Department of Mathematics, Celestijnenlaan
200B box 2400,
${}^{{\dagger}^{4}}$ BE-3001 Leuven, Belgium<EMAIL_ADDRESS>
Received January 05, 2016, in final form June 01, 2016; Published online June
13, 2016

In this paper we present multidimensional analogues of both the continuous-
and discrete-time Toda lattices. The integrable systems that we consider here
have two or more space coordinates. To construct the systems, we generalize
the orthogonal polynomial approach for the continuous and discrete Toda
lattices to the case of multiple orthogonal polynomials.

multiple orthogonal polynomials; orthogonal polynomials; recurrence relations;
Toda equation; discrete integrable system; Toda lattice

42C05; 37K10; 39A14; 65Q10

##  1 Introduction

The continuous-time Toda lattice [53, 54]

$\displaystyle\dot{a}_{n}(t)=a_{n}(t)(b_{n-1}(t)-b_{n}(t)),$
$\displaystyle\dot{b}_{n}(t)=a_{n}(t)-a_{n-1}(t),\qquad a_{n}>0,\qquad
t\in\mathbb{R}_{+},\qquad n\in{\mathbb{Z}}_{+},$ (1.1)

and the discrete-time Toda lattice [31, 50]

$\displaystyle A_{n}^{t+1}+B_{n}^{t+1}=A_{n}^{t}+B_{n+1}^{t},$ $\displaystyle
A_{n-1}^{t+1}B_{n}^{t+1}=A_{n}^{t}B_{n}^{t},\qquad t,n\in{\mathbb{Z}}_{+},$
(1.2)

have appeared in physical and mathematical models quite a while ago and are
still attracting the interest of many researchers in the field. For instance,
generalizations of the Toda lattice have been considered from various points
of view (e.g., [33, 48]). One of the interesting aspects of the continuous-
and discrete-time Toda lattice is that orthogonal polynomials (OPs) appear as
eigenfunctions of their Lax pairs [11, 50], which means that spectral
transformations of OPs describe the flow of the Toda lattice. In addition,
several integrable systems have been shown to be related to descendants of OPs
through their spectral transformations [10, 22, 23, 24]. Applying this
spectral transformation technique to a new class of OPs, novel integrable
systems have been exploited [1, 50, 51]. Recall (see [36, 41]) that the system
(1.1) is managed by the evolution

$\displaystyle d\mu(x,t)=e^{-xt}d\mu(x),$ (1.3)

and the functions $a_{n}(t)$ and $b_{n}(t)$ appear as the coefficients of the
three-term recurrence relation

$\displaystyle
xP_{n}(x;t)=P_{n+1}(x;t)+b_{{n}}(t)P_{{n}}(x;t)+a_{{n}}(t)P_{{n}-1}(x;t),\qquad
n\in{\mathbb{Z}}_{+},$ (1.4)

where $P_{{n}}(x;t)$ are monic polynomials in the variable $x$, orthogonal
with respect to $d\mu(x,t)$. Thus, the direct and inverse spectral
transformations $\\{a_{{n}}(t),b_{{n}}(t)\\}\rightleftarrows d\mu(t,x)$ along
with the evolution (1.3) solve the Cauchy problem for (1.1).

In this paper, motivated by these results, we aim to explore a new
generalization of the Toda lattice by developing the spectral transformations
of multiple orthogonal polynomials (m-OPs) [6, 32], which originated from the
theory of Hermite–Padé approximants [44], which were introduced by Hermite
[30] in connection with his proof of the transcendence of $e$. On top of that,
m-OPs were recently found to have applications in many areas such as random
matrices [9, 19] and the theory of difference operators on lattices [7, 13].

A lattice of multiple orthogonal polynomials $P_{\vec{n}}$,
$\vec{n}\in{\mathbb{Z}}_{+}^{r}$ is defined as a set of polynomials of degree
$|\vec{n}|=n_{1}+\cdots+n_{r}$, satisfying the orthogonality relations

$\displaystyle\mathcal{L}_{j}[x^{i}P_{\vec{n}}(x)]=\int
x^{i}P_{\vec{n}}(x)d\mu_{j}(x)=0,\qquad i=0,\dots,n_{j}-1,\qquad 1\leq j\leq
r.$ (1.5)

This definition gives $|\vec{n}|$ homogeneous and linear equations for the
$|\vec{n}|+1$ coefficients of the polynomial $P_{\vec{n}}$. In fact, we can
choose $P_{\vec{n}}$ to be monic and, so, it leads to a linear system of
$|\vec{n}|$ unknown coefficients, which has a unique solution if and only if
the corresponding determinant is not vanishing. Evidently, it is not always
the case that the determinant is nonzero. Therefore, there is no guarantee
that for a given multi-index one can find the corresponding multiple
orthogonal polynomial. In the case of uniqueness the multi-index $\vec{n}$ is
called normal. Hence, for normal indices the monic polynomial $P_{\vec{n}}$
can be determined. If all multi-indices $\vec{n}\in{\mathbb{Z}}_{+}^{r}$ on
the lattice are normal then the system of measures $\\{\mu_{j}\\}_{j=1}^{r}$
(or functionals $\\{\mathcal{L}_{j}\\}_{j=1}^{r}$) generating the lattice of
polynomials $\\{P_{\vec{n}}\\}$ in (1.5) is called a perfect system. In other
words, for a perfect system the polynomial $P_{\vec{n}}$ is well defined for
any $\vec{n}\in{\mathbb{Z}}_{+}^{r}$. This notion has been introduced by K.
Mahler [37] in relation to the diophantine approximation.

It is noteworthy that m-OPs become ordinary OPs if we take $r=1$. Moreover,
one of the properties inherited from ordinary OPs is that m-OPs satisfy the
following nearest-neighbor recurrence relations [32, 55], which actually
generalizes the three-term recurrence relation for OPs (1.4):

$\displaystyle
xP_{\vec{n}}(x)=P_{\vec{n}+\vec{e}_{i}}(x)+b_{\vec{n},i}P_{\vec{n}}(x)+\sum_{k=1}^{r}a_{\vec{n},k}P_{\vec{n}-\vec{e}_{k}}(x),\qquad
i=1,\dots,r,\qquad\vec{n}\in{\mathbb{Z}}_{+}^{r},$ (1.6)

where $\vec{e}_{i}=(0,\dots,0,1,0,\dots,0)$ is the $i$-th vector of the
standard basis in ${\mathbb{Z}}_{+}^{r}$. Unlike in the ordinary case, the
coefficients of the recurrence relations $\\{a_{\vec{n},i},b_{\vec{n},i}\\}$
for m-OPs are not independent but are required to satisfy the difference
equations [55]

$\displaystyle\sum_{k=1}^{r}(a_{\vec{n}+\vec{e}_{j},k}-a_{\vec{n}+\vec{e}_{i},k})=(b_{\vec{n},j}-b_{\vec{n},i})(b_{\vec{n}+\vec{e}_{j},i}-b_{\vec{n},i}),$
$\displaystyle\frac{a_{\vec{n},i}}{a_{\vec{n}+\vec{e}_{j},i}}=\frac{b_{\vec{n}-\vec{e}_{i},j}-b_{\vec{n}-\vec{e}_{i},i}}{b_{\vec{n},j}-b_{\vec{n},i}},\qquad
1\leq i\neq j\leq r.$ (1.7)

Observe that (1.7) is equivalent to the equations in [55, Theorem 3.2]:
interchanging $i$ and $j$ in the first equation of (1.7) gives (3.6) from [55]
and the determinant in (3.7) from [55] is the same as the right hand side of
the first equation in our (1.7).

Finally, we are in a position to state the main results of the paper.

######  Theorem 1.1.

Suppose that the system of measures

$\displaystyle d\mu_{j}(x,t)=e^{-tx}d\mu_{j}(x,0),\qquad 1\leq j\leq r,$ (1.8)

generates m-OPs with normal indices $\vec{n}$ and
$\\{\vec{n}\pm\vec{e}_{k}\\}_{k=1}^{r}$ in a neighborhood of $t=0$. Then the
recurrence coefficients from (1.6) satisfy locally the equations

$\displaystyle\dot{a}_{\vec{n},k}=a_{\vec{n},k}[b_{\vec{n}-\vec{e}_{k},k}-b_{\vec{n},k}],$
$\displaystyle\dot{b}_{\vec{n},k}=\sum_{j=1}^{r}(a_{\vec{n},j}-a_{\vec{n}+\vec{e}_{k},j}),\qquad
1\leq k\leq r.$ (1.9)

If the system $\\{d\mu_{j}(x,t)\\}_{j=1}^{r}$ in (1.8) is perfect for
$t\geqslant 0$, then the Cauchy problem for (1.9) has a global solution which
can be obtained by the direct and inverse spectral transformations
$\\{a_{\vec{n},k}(t),b_{\vec{n},k}(t)\\}\rightleftarrows\\{d\mu_{k}(x,t)\\}$
and the evolution (1.8).

######  Remark 1.2.

We would like to emphasize here that the system (1.9) is solvable as long as
we choose the initial values subject to (1.7), which in turn means that we
have to start with a perfect system of measures. Once again, such systems
generate coefficients that satisfy (1.7). Moreover, for the global solution
the Toda dynamics preserve the stationary equations (1.7). Therefore, (1.9)
forms a “weakly integrable” system. Similar phenomena occur for the quantum
Hamiltonians associated with classical multiple orthogonal polynomials [39,
40, 42]. Also, it should be noted that there are other multidimensional
integrable systems studied in the literature. For instance, see [38]; see also
[8] and references therein.

######  Remark 1.3.

There are two well-known perfect systems of measures. One of them is called an
Angelesco system [5] and is formed by measures with supports on disjoint
intervals. The perfectness of an Angelesco system immediately follows from
(1.5) and the properties of zeros of OPs. The other one is called a Nikishin
system [43]. The measures from a Nikishin system have the same support but
some extra conditions (analytic properties of the weight functions) need to
hold. Details of the definition of Nikishin systems and the proof of their
perfectness can be found in [26].

######  Remark 1.4.

The system $\\{d\mu_{j}\\}_{j=1}^{r}$ consists of the spectral measures of the
$r$ marginal one-dimensional difference operators defined by the three-term
recurrence relations (1.4) with coefficients $a_{n}^{(j)}:=a_{ne_{j},j}$,
$b_{n}^{(j)}:=b_{ne_{j},j}$, $j=1,\dots,r$. These marginal spectral measures
can be taken as spectral data for the multidimensional difference operator
defined by the recurrence relations (1.6) with the coefficients
$\\{a_{\vec{n},k},b_{\vec{n},k}\\}$. Therefore, the direct and inverse
spectral transformations
$\\{a_{\vec{n},k}(t),b_{\vec{n},k}(t)\\}\rightleftarrows\\{d\mu_{k}(t,x)\\}$
reduce to the well-known direct and inverse spectral problems for OPs together
with a scheme to solve a boundary value problem (BVP) for the discrete
integrable system (1.7); see [14, 27] for more details on this matter.

The proof of Theorem 1.1 and properties of m-OPs when the corresponding
measures evolve as in (1.8) are presented in Section 2 (see Subsections
2.1–2.3).

It turns out that an important place on the lattice of m-OPs (1.5)–(1.6) is
the diagonal. Namely, the diagonal sequence $\\{q_{N}\\}$ of multiple
orthogonal polynomials [45] (also called $r$-orthogonal polynomials [16, 17,
25]) that is generated in the following manner

$\displaystyle q_{N}=P_{(n,\dots,n)+\sum\limits_{j=1}^{k}e_{j}},\qquad
N=:nr+k,\qquad k=0,\dots,r-1.$ (1.10)

The diagonal m-OPs appear as the common denominator of the convergents of the
underlying Jacobi–Perron vector continued fraction [18, 44], which is a
functional analog of the continued fraction introduced by Jacobi [34, 35] and
Perron [47] in their approach to find a characterization of irrationals of
orders higher than quadratic. As a matter of fact, the polynomial sequence
$\\{q_{N}\\}$ satisfies the step-line recurrence relation [6, 44]

$\displaystyle
xq_{N}(x)=q_{N+1}(x)+\beta_{N}q_{N}(x)+\sum_{k=1}^{r}\alpha_{N}^{(k)}q_{N-k}(x),\qquad
n\in{\mathbb{Z}}_{+}.$ (1.11)

If the system of measures $\\{d\mu_{j}\\}_{j=1}^{r}$ evolves subject to (1.8),
then (see [2]) the equations for the coefficients from (1.11) can be written
in the Lax pair form

$\displaystyle\dot{L}=[L,(L)_{-}],\qquad
L:=L\big{(}1,\beta_{N},\alpha_{N}^{(1)},\dots,\alpha_{N}^{(r)}\big{)},$ (1.12)

where $(X)_{-}$ denotes the strictly lower part of the matrix $X$ and $L$ is a
$(r+2)$-banded (lower Hessenberg) semi-infinite matrix with 1’s on the upper
diagonal, $\\{\beta_{N}\\}$ on the main diagonal and
$\\{\alpha_{N}^{(k)}\\}_{k=1}^{r}$ building the lower diagonals. In Subsection
2.4 we show relations between two generalizations of the Toda equations, which
are basically the relations between (1.9) and (1.12).

In the same spirit as it is done for the continuous-time Toda equation, we
obtain in Section 3 a discrete analogue of (1.8) and (1.9).

######  Theorem 1.5.

Suppose that the system of measures

$\displaystyle d\mu_{j}(x,t)=x^{t}d\mu_{j}(x,0),\qquad
t\in{\mathbb{Z}}_{+},\qquad 1\leq j\leq r,$

generates the m-OPs $\\{P_{\vec{n}}^{t}(x)\\}$ with normal indices $\vec{n}$
and $\\{\vec{n}\pm\vec{e}_{k}\\}_{k=1}^{r}$. Then the following nonlinear
difference system on a semi-infinite lattice

$\displaystyle
A_{\vec{n},j}^{t+1}+\sum_{k=1}^{r}B_{\vec{n},k}^{t+1}=A_{\vec{n},j}^{t}+\sum_{k=1}^{r}B_{\vec{n}+\vec{e_{j}},k}^{t},$
$\displaystyle
A_{\vec{n}-\vec{e}_{j},j}^{t+1}B_{\vec{n},j}^{t+1}=A_{\vec{n},j}^{t}B_{\vec{n},j}^{t},\qquad
1\leq j\leq r,$ (1.13)

can be solved locally with respect to the space variable. Moreover, if the
system of measures is perfect for any discrete time $t\in{\mathbb{Z}}_{+}$,
then the solution exists globally and it is given by the formulas

$\displaystyle
A_{\vec{n},j}^{t}=-\frac{P_{\vec{n}+\vec{e}_{j}}^{t}(0)}{P_{\vec{n}}^{t}(0)},\qquad
B_{\vec{n},k}^{t}=-a_{\vec{n},k}^{t}\frac{P_{\vec{n}-\vec{e}_{k}}^{t}(0)}{P_{\vec{n}}^{t}(0)}.$

To get to (1.13) we present two approaches. The first one is based on
Christoffel and Geronimus transformations, which are known to be discrete
analogues of Darboux transformations, see Subsection 3.2. The second method is
obtained by following the consistency approach from [20] and [49]. In
particular, the Lax pair we get for the second method is a certain adaptation
of the one from [49] to our setting. The latter approach allows us to present
the discrete-time Toda equations (1.13) and the consistency equations (1.7) in
a unified fashion in the form of Lax pairs commutation relations, see
Subsection 3.4. It is worth mentioning that both the approaches are related to
a generalization of the quotient-difference algorithm for the Padé table.
Moreover, our approach complements some earlier attempts to develop the q-d
algorithm [56] related to multiple orthogonal polynomials and puts it into the
context of discrete integrable systems.

Also, in Subsection 3.3 we connect the discrete-time multidimensional Toda
equations (1.13) and the discrete time analogue of Toda chain equation (1.12)
for the recurrence coefficients of the diagonal sequence of m-OPs (1.11).
Using the formulas for multiple Laguerre polynomials [12, 55] we give the
initial data for the explicit solution of the multidimensional Toda lattice in
both cases. For the continuous time it can be found in Subsection 2.5 and for
the discrete time in Subsection 3.5.

##  2 The continuous-time higher analogues of the Toda lattice

###  2.1 Preliminaries

Let us briefly go over the basic definitions, see [32, Chapter 23] for
details. The m-OPs $P_{\vec{n}}$ defined in (1.5) are called multiple
orthogonal (or Hermite–Padé) polynomials of type II. In what follows, we
suppose that $P_{\vec{n}}$ are monic so that the m-OPs are uniquely
determined. It is useful to consider the dual construction. More precisely,
type I multiple orthogonal polynomials $(C_{\vec{n},1},\dots,C_{\vec{n},r})$
are such that $C_{\vec{n},j}$ is a polynomial of degree at most ${n}_{j}-1$
with the orthogonality relations

$\displaystyle\sum^{r}_{j=1}\int x^{k}C_{\vec{n},j}(x)d\mu_{j}(x)=0,\qquad
k=0,1,\dots,|\vec{n}|-2,$

and the normalization

$\displaystyle\sum^{r}_{j=1}\int
x^{|\vec{n}|-1}C_{\vec{n},j}(x)d\mu_{j}(x)=1,\qquad
C_{\vec{n},j}(x)=\kappa_{\vec{n},j}x^{n_{j}-1}+\cdots.$

We assume that the $r$ measures $\mu_{1},\dots,\mu_{r}$ are all absolutely
continuous with respect to a measure $\mu$ and that
$d\mu_{j}(x)=w_{j}(x)d\mu(x)$ and we will use the following notation

$\displaystyle Q_{\vec{n}}(x)=\sum^{r}_{j=1}C_{\vec{n},j}(x)w_{j}(x).$

It is easy to check that the type I and type II multiple orthogonal
polynomials form a biorthonormal system in the sense that the functions
$Q_{\vec{n}}$ generated by type I multiple orthogonal polynomials satisfy

$\displaystyle\int P_{\vec{n}}(x)Q_{\vec{m}}(x)d\mu(x)=\begin{cases}0&\text{if
}|\vec{m}|\leq|\vec{n}|,\\\ 0&\text{if }|\vec{n}|\leq|\vec{m}|-2,\\\
1&\text{if }\vec{m}=\vec{n}+\vec{e}_{k}\text{ for }1\leq k\leq r.\end{cases}$
(2.1)

Furthermore, for the type II polynomials we have (1.6):

$\displaystyle
xP_{\vec{n}}(x)=P_{\vec{n}+\vec{e}_{k}}(x)+b_{\vec{n},k}P_{\vec{n}}(x)+\sum^{r}_{j=1}a_{\vec{n},j}P_{\vec{n}-\vec{e}_{j}}(x),$
(2.2)

and using the type I polynomials we get

$\displaystyle b_{\vec{n},k}=\int
xP_{\vec{n}}Q_{\vec{n}+\vec{e}_{k}}(x)d\mu(x),\qquad a_{\vec{n},j}=\frac{\int
x^{n_{j}}P_{\vec{n}}(x)d\mu_{j}(x)}{\int
x^{n_{j}-1}P_{\vec{n}-\vec{e}_{k}}(x)d\mu_{j}(x)}.$ (2.3)

For type I we have a similar recurrence relation

$\displaystyle
xQ_{\vec{n}}(x)=Q_{\vec{n}-\vec{e}_{k}}(x)+b_{\vec{n}-\vec{e}_{k},k}Q_{\vec{n}}(x)+\sum^{r}_{j=1}a_{\vec{n},j}Q_{\vec{n}+\vec{e}_{j}}(x).$
(2.4)

We need to point out two useful relations here. If we multiply (2.2) by
$Q_{\vec{n}}(x)$ and integrate, then the biorthogonality gives

$\displaystyle\int
xP_{\vec{n}}(x)Q_{\vec{n}}(x)d\mu(x)=\sum_{j=1}^{r}a_{\vec{n},j}.$ (2.5)

The orthogonality properties of $P_{\vec{n}}$ imply that

$\displaystyle\int P_{\vec{n}}(x)Q_{\vec{n}+\vec{e}_{k}}(x)d\mu(x)=\int
P_{\vec{n}}C_{\vec{n}+\vec{e}_{k},k}(x)d\mu_{k}(x),$

so we arrive at

$\displaystyle 1=\kappa_{\vec{n}+\vec{e}_{k},k}\int
P_{\vec{n}}(x)x^{n_{k}}d\mu_{k}(x).$ (2.6)

###  2.2 Time dynamics

In this subsection we present a proof of Theorem 1.1. To this end, we consider
a perfect system of measures $(\mu_{1}(x,t),\dots,\mu_{r}(x,t))$ that depend
on time as in (1.8), i.e., $d\mu_{k}(x,t):=e^{-tx}d\mu_{k}(x)$ for $1\leq
k\leq r$. Now the corresponding type I and type II multiple orthogonal
polynomials depend on time as well, that is, we have $Q_{\vec{n}}(x;t)$ and
$P_{\vec{n}}(x;t)$. Our goal is to show that the relations (2.3) evaluated at
the moment $t$

$\displaystyle b_{\vec{n},k}(t)=\int
xP_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x),\qquad
a_{\vec{n},k}(t)=\frac{\int x^{n_{k}}P_{\vec{n}}(x;t)e^{-xt}d\mu_{k}(x)}{\int
x^{n_{k}-1}P_{\vec{n}-\vec{e}_{k}}(x)e^{-xt}d\mu_{k}(x)},$

imply multiple Toda (m-Toda) equations (1.9), i.e., for $1\leq k\leq r$

$\displaystyle\dot{a}_{\vec{n},k}=a_{\vec{n},k}[b_{\vec{n}-\vec{e}_{k},k}-b_{\vec{n},k}],$
(2.7)

and

$\displaystyle\dot{b}_{\vec{n},k}=\sum_{j=1}^{r}(a_{\vec{n},j}-a_{\vec{n}+\vec{e}_{k},j}).$
(2.8)

From the biorthogonality (2.1) we get

$\displaystyle\int
P_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=1.$

If we take derivatives with respect to $t$, then this gives

$\displaystyle\int\dot{P}_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)+\int
P_{\vec{n}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)$
$\displaystyle\qquad{}=\int
xP_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x).$

Observe that $\dot{P}_{\vec{n}}$ is a polynomial of degree at most
$|\vec{n}|-1$ since ${P}_{\vec{n}}$ is a monic polynomial, hence the
orthogonality for type I multiple orthogonal polynomials gives

$\displaystyle\int\dot{P}_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=0.$

The orthogonality for type II multiple orthogonal polynomials gives

$\displaystyle\int
P_{\vec{n}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=\int
P_{\vec{n}}(x;t)\dot{A}_{\vec{n}+\vec{e}_{k},k}(x;t)e^{-xt}d\mu_{k}(x).$ (2.9)

If we use
$C_{\vec{n}+\vec{e}_{k},k}(x;t)=\kappa_{\vec{n}+\vec{e}_{k},k}(t)x^{n_{k}}+\cdots$,
then we find

$\displaystyle\int
P_{\vec{n}}(x;t)\dot{A}_{\vec{n}+\vec{e}_{k},k}(x;t)e^{-xt}d\mu_{k}(x)=\kappa_{\vec{n}+\vec{e}_{k},k}\int
x^{n_{k}}P_{\vec{n}}(x;t)e^{-xt}d\mu_{k}(x)=\frac{\dot{\kappa}_{\vec{n},j}}{\kappa_{\vec{n}+\vec{e}_{k},j}},$

where the last equality follows from (2.5). Combining these results with (2.3)
already gives

$\displaystyle
b_{\vec{n},k}(t)=\frac{\dot{\kappa}_{\vec{n}+\vec{e}_{k},k}}{\kappa_{\vec{n}+\vec{e}_{k},k}}.$
(2.10)

Taking derivatives in (2.6) gives

$\displaystyle\frac{\dot{\kappa}_{\vec{n},k}}{\kappa_{\vec{n},k}}-\frac{\dot{\kappa}_{\vec{n}+\vec{e}_{k},k}}{\kappa_{\vec{n}+\vec{e}_{k},k}}=\frac{\dot{a}_{\vec{n},k}}{a_{\vec{n},k}},$

which implies

$\displaystyle
b_{\vec{n}-\vec{e}_{k},k}(t)-b_{\vec{n},k}(t)=\frac{\dot{a}_{\vec{n}+\vec{e}_{k},k}}{a_{\vec{n}+\vec{e}_{k},k}},$

giving (2.7). Thus the first relation in (1.9) is proved.

For the second relation in (1.9) we use (2.3) to find

$\displaystyle b_{\vec{n},k}(t)=\int
xP_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x).$

Taking the derivative with respect to $t$ gives

$\displaystyle\dot{b}_{\vec{n},k}(t)=\int
x\dot{P}_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)$
$\displaystyle\hphantom{\dot{b}_{\vec{n},k}(t)=}{}+\int
xP_{\vec{n}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)-\int
x^{2}P_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x).$

If we use (2.4) then we find

$\displaystyle\int
x\dot{P}_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)$
$\displaystyle\qquad{}=\int\dot{P}_{\vec{n}}(x;t)\left(Q_{\vec{n}}(x;t)+b_{\vec{n},k}Q_{\vec{n}+\vec{e}_{k}}(x;t)+\sum^{r}_{j=1}a_{\vec{n}+\vec{e}_{k},j}Q_{\vec{n}+\vec{e}_{j}+\vec{e}_{k}}(x;t)\right)e^{-xt}d\mu(x).$

Since $\dot{P}_{\vec{n}}$ is of degree at most $|\vec{n}|-1$, the
orthogonality for type I multiple orthogonal polynomials gives

$\displaystyle\int
x\dot{P}_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=\int\dot{P}_{\vec{n}}(x;t)Q_{\vec{n}}(x;t)e^{-xt}d\mu(x).$

Taking the derivative of

$\displaystyle\int{P}_{\vec{n}}(x;t)Q_{\vec{n}}(x;t)e^{-xt}d\mu(x)=0$

leads to

$\displaystyle\int\dot{P}_{\vec{n}}(x;t)Q_{\vec{n}}(x;t)e^{-xt}d\mu(x)=\int
x{P}_{\vec{n}}(x;t)Q_{\vec{n}}(x;t)e^{-xt}d\mu(x)=\sum^{r}_{j=1}a_{\vec{n},j}(t),$

hence

$\displaystyle\int
x\dot{P}_{\vec{n}}(x;t)Q_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=\sum^{r}_{j=1}a_{\vec{n},j}(t).$

If we use (2.2) then we find

$\displaystyle\int
x{P}_{\vec{n}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)$
$\displaystyle\qquad{}=\int\left(P_{\vec{n}+\vec{e}_{k}}(x;t)+b_{\vec{n},k}(t)P_{\vec{n}}(x;t)+\sum^{r}_{j=1}a_{\vec{n},j}P_{\vec{n}-\vec{e}_{j}}(x;t)\right)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x).$

The orthogonality of type II multiple orthogonal polynomials gives

$\displaystyle\int
x{P}_{\vec{n}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)$
$\displaystyle\qquad{}=b_{\vec{n},k}\int{P}_{\vec{n}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)+\sum^{r}_{j=1}a_{\vec{n},j}\int{P}_{\vec{n}-\vec{e}_{j}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x).$

From (2.9) and (2.10) we recall that

$\displaystyle\int{P}_{\vec{n}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=b_{\vec{n},k}.$

If we take the derivative of

$\displaystyle\int{P}_{\vec{n}-\vec{e}_{j}}(x;t){Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=0,$

then we find

$\displaystyle\int{P}_{\vec{n}-\vec{e}_{j}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=\int
x{P}_{\vec{n}-\vec{e}_{j}}(x;t){Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x),$

and after using (2.2) and the biorthogonality, we find

$\displaystyle\int{P}_{\vec{n}-\vec{e}_{j}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=1.$

This gives

$\displaystyle\int
x{P}_{\vec{n}}(x;t)\dot{Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)={b}^{2}_{\vec{n},k}+\sum^{r}_{j=1}a_{\vec{n},j}.$

Finally use (2.2), (2.4) and the biorthogonality to find

$\displaystyle\int
x^{2}{P}_{\vec{n}}(x;t){Q}_{\vec{n}+\vec{e}_{k}}(x;t)e^{-xt}d\mu(x)=\sum^{r}_{j=1}a_{\vec{n}+\vec{e}_{k},j}+{b}^{2}_{\vec{n},k}+\sum^{r}_{j=1}a_{\vec{n},j}.$

Combining all these results then gives

$\displaystyle\dot{b}_{\vec{n},k}=\sum^{r}_{j=1}a_{\vec{n},j}-\sum^{r}_{j=1}a_{\vec{n}+\vec{e}_{k},j},$

which is the same as (2.8). This finishes the proof of Theorem 1.1.

###  2.3 Time dependent m-OPs

Introducing the moments defined by the functional (1.5)

$\displaystyle\mu_{i,j}:=\mathcal{L}_{j}[x^{i}],\qquad j=1,\dots,r,$ (2.11)

we have a determinant expression for m-OPs

$\displaystyle
P_{\vec{n}}(x)=\frac{1}{\tau_{\vec{n}}}\begin{vmatrix}\mu_{0,1}&\cdots&\mu_{n_{1}-1,1}&\cdots&\mu_{0,r}&\cdots&\mu_{n_{r}-1,r}&1\\\
\mu_{1,1}&\cdots&\mu_{n_{1},1}&\cdots&\mu_{1,r}&\cdots&\mu_{n_{r},r}&x\\\
\vdots&\cdots&\vdots&\vdots&\cdots&\vdots&\cdots&\vdots\\\
\mu_{|\vec{n}|,1}&\cdots&\mu_{|\vec{n}|+n_{1}-1,1}&\cdots&\mu_{|\vec{n}|,r}&\cdots&\mu_{|\vec{n}|+n_{r}-1,r}&x^{|\vec{n}|}\end{vmatrix},$
(2.12)

with

$\displaystyle\tau_{\vec{n}}=\begin{vmatrix}\mu_{0,1}&\cdots&\mu_{n_{1}-1,1}&\cdots&\mu_{0,r}&\cdots&\mu_{n_{r}-1,r}\\\
\mu_{1,1}&\cdots&\mu_{n_{1},1}&\cdots&\mu_{1,r}&\cdots&\mu_{n_{r},r}\\\
\vdots&\cdots&\vdots&\cdots&\vdots&\cdots&\vdots\\\
\mu_{|\vec{n}|-1,1}&\cdots&\mu_{|\vec{n}|+n_{1}-2,1}&\cdots&\mu_{|\vec{n}|-1,r}&\cdots&\mu_{|\vec{n}|+n_{r}-2,r}\end{vmatrix}.$

From this expression, it is easy to find that the multi-index
$\vec{n}\in\mathbb{Z}_{+}^{n}$ is normal iff $\tau_{\vec{n}}\neq 0$ and this
is assumed to hold in what follows. Using (2.12), (1.6) and the orthogonality
relation (1.5), the following determinant expression of the recurrence
coefficients $\\{a_{\vec{n},j},b_{\vec{n},j}\\}$ is directly verified

$\displaystyle
b_{\vec{n},j}=\frac{\sigma_{\vec{n}+\vec{e}_{j}}}{\tau_{\vec{n}+\vec{e}_{j}}}-\frac{\sigma_{\vec{n}}}{\tau_{\vec{n}}},\qquad
a_{\vec{n},j}=\frac{\tau_{\vec{n}+\vec{e}_{j}}\tau_{\vec{n}-\vec{e}_{j}}}{\tau_{\vec{n}}^{2}},$
(2.13)

with

$\displaystyle\sigma_{\vec{n}}=\begin{vmatrix}\mu_{0,1}&\cdots&\mu_{n_{1}-1,1}&\cdots&\mu_{0,r}&\cdots&\mu_{n_{r}-1,r}\\\
\vdots&\cdots&\vdots&\cdots&\vdots&\cdots&\vdots\\\
\mu_{|\vec{n}|-2,1}&\cdots&\mu_{|\vec{n}|+n_{1}-3,1}&\cdots&\mu_{|\vec{n}|-2,r}&\cdots&\mu_{|\vec{n}|+n_{r}-3,r}\\\
\mu_{|\vec{n}|,1}&\cdots&\mu_{|\vec{n}|+n_{1}-1,1}&\cdots&\mu_{|\vec{n}|,r}&\cdots&\mu_{|\vec{n}|+n_{r}-1,r}\end{vmatrix}.$

As was already mentioned, the spectral transformation plays a central role in
finding the corresponding integrable systems. In order to mimic the classical
scheme, we will consider the spectral transformation of m-OPs first. Let us
introduce the 1-parameter deformation of the moments (2.11) as follows

$\displaystyle\frac{d}{dt}\mu_{i,j}=-\mu_{i+1,j},\qquad j=1,\dots,r.$ (2.14)

This transformation can also be interpreted in terms of the linear functionals
(1.5):

$\displaystyle\frac{d}{dt}\mathcal{L}_{j}[\cdot]=-\mathcal{L}_{j}[x\cdot],\qquad
j=1,\dots,r.$ (2.15)

Notice that the 1-parameter deformation (2.14) (or (2.15)) coincides with that
of ordinary OPs in the case $r=1$. Using the determinant expression of m-OPs
(2.12), the spectral transformation of m-OPs can be constructed.

######  Theorem 2.1.

If the linear functionals with the condition (2.15) $($or equivalently the
moments $\\{\mu_{i,j}\\}$ with (2.14)$)$ are given, then the following
relation for the corresponding m-OPs holds

$\displaystyle\frac{d}{dt}P_{\vec{n}}(x)=\sum_{k=1}^{r}a_{\vec{n},k}P_{\vec{n}-\vec{e}_{k}}(x),$
(2.16)

where $\\{a_{\vec{n},j}\\}$ are the coefficients of the recurrence relation in
(1.6).

###### Proof.

To begin with, introduce the notation

$\displaystyle\tau_{\vec{n},x}:=\big{|}0^{n_{1}},1^{n_{1}},\dots,(n_{1}-1)^{n_{1}},\dots,0^{n_{r}},\dots,(n_{r}-1)^{n_{r}},x\big{|}$
$\displaystyle\hphantom{\tau_{\vec{n},x}}{}=\begin{vmatrix}\mu_{0,1}&\mu_{1,1}&\cdots&\mu_{n_{1}-1,1}&\cdots&\mu_{0,r}&\cdots&\mu_{n_{r}-1,r}&1\\\
\mu_{1,1}&\mu_{2,1}&\cdots&\mu_{n_{1}+1,1}&\cdots&\mu_{1,r}&\cdots&\mu_{n_{r},r}&x\\\
\vdots&\vdots&\cdots&\vdots&\cdots&\vdots&\cdots&\vdots&\vdots\\\
\mu_{|\vec{n}|,1}&\mu_{|\vec{n}|+1,1}&\cdots&\mu_{|\vec{n}|+n_{1}-1,1}&\cdots&\mu_{|\vec{n}|,r}&\cdots&\mu_{|\vec{n}|+n_{r}-1,r}&x^{|\vec{n}|}\end{vmatrix}.$

Then $\tau_{\vec{n}}$ and $P_{\vec{n}}(x)$ can be rewritten as

$\displaystyle
P_{\vec{n}}(x)=\frac{\tau_{\vec{n},x}}{\tau_{\vec{n}}},\qquad\tau_{\vec{n}}=\big{|}0^{n_{1}},\dots,(n_{1}-1)^{n_{1}},\dots,0^{n_{r}},\dots,(n_{r}-1)^{n_{r}-1},e\big{|},$

where $e=(0,\dots,0,1)^{T}$. From the relation (2.14), it is easy to find that
$\tau_{\vec{n}}$ and $\tau_{\vec{n},x}$ are Wronskian matrices, which amounts
to

$\displaystyle-\frac{d}{dt}\tau_{\vec{n}}=\sum_{k=1}^{r}\big{|}0^{n_{1}},\dots,(n_{k}-2)^{n_{k}},n_{k}^{n_{k}},0^{n_{k+1}},\dots,(n_{r}-1)^{n_{r}-1},e\big{|},$
$\displaystyle-\frac{d}{dt}\tau_{\vec{n},x}=\sum_{k=1}^{r}\big{|}0^{n_{1}},\dots,(n_{k}-2)^{n_{k}},n_{k}^{n_{k}},0^{n_{k+1}},\dots,(n_{r}-1)^{n_{r}-1},x\big{|}.$

Using these notations and relations, we can calculate the derivative of
$P_{\vec{n}}(x)$:

$\displaystyle-\frac{d}{dt}P_{\vec{n}}(x)=-\frac{d}{dt}\frac{\tau_{\vec{n},x}}{\tau_{\vec{n}}}=-\frac{\frac{d}{dt}\tau_{\vec{n},x}\tau_{\vec{n}}-\tau_{\vec{n},x}\frac{d}{dt}\tau_{\vec{n}}}{\tau_{\vec{n}}^{2}}$
$\displaystyle\hphantom{-\frac{d}{dt}P_{\vec{n}}(x)}{}=-\frac{1}{\tau_{\vec{n}}^{2}}\sum_{k=1}^{r}\Bigl{(}\big{|}\dots,(n_{k}-2)^{n_{k}},n_{k}^{n_{k}},\dots,x\big{|}\big{|}\dots,(n_{k}-2)^{n_{k}},(n_{k}-1)^{n_{k}},\dots,e\big{|}$
$\displaystyle\hphantom{-\frac{d}{dt}P_{\vec{n}}(x)=}{}-\big{|}\dots,(n_{k}-2)^{n_{k}},(n_{k}-1)^{n_{k}},\dots,x\big{|}\big{|}\dots,(n_{k}-2)^{n_{k}},n_{k}^{n_{k}},\dots,e\big{|}\Bigr{)}$
$\displaystyle\hphantom{-\frac{d}{dt}P_{\vec{n}}(x)}{}=\frac{1}{\tau_{\vec{n}}^{2}}\sum_{k=1}^{r}\big{|}\dots,(n_{k}-2)^{n_{k}},\dots,x,e\big{|}\big{|}\dots,(n_{k}-2)^{n_{k}},(n_{k}-1)^{n_{k}},n_{k}^{n_{k}},\dots\big{|}$
$\displaystyle\hphantom{-\frac{d}{dt}P_{\vec{n}}(x)}{}=\sum_{k=1}^{r}\frac{\tau_{\vec{n}-\vec{e}_{k},x}\tau_{\vec{n}+\vec{e}_{k}}}{\tau_{\vec{n}}^{2}}=\sum_{k=1}^{r}\frac{\tau_{\vec{n}-\vec{e}_{k}}\tau_{\vec{n}+\vec{e}_{k}}}{\tau_{\vec{n}}^{2}}\frac{\tau_{\vec{n}-\vec{e}_{k},x}}{\tau_{\vec{n}-\vec{e}_{k}}}.$
(2.17)

In the calculation of (2.17), we have used the Plücker relation, a well-known
identity for determinants

$\displaystyle|\dots,a,b||\dots,c,d|-|\dots,a,c||\dots,b,d|+|\dots,a,d||\dots,b,c|=0,$

where $a$, $b$, $c$, $d$ are arbitrary column vectors of appropriate size.
Finally, comparing the result with (2.13), we arrive at (2.16). This completes
the proof. ∎

We thus have obtained $r+1$ linear equations where m-OPs appear as an
eigenfunction (we shall refer to this as a “Lax set”)

$\displaystyle
xP_{\vec{n}}(x)=P_{\vec{n}+\vec{e}_{j}}(x)+b_{\vec{n},j}P_{\vec{n}}(x)+\sum_{k=1}^{r}a_{\vec{n},k}P_{\vec{n}-\vec{e}_{k}}(x),\qquad
j=1,\dots,r,$
$\displaystyle\frac{d}{dt}P_{\vec{n}}(x)=\sum_{k=1}^{r}a_{\vec{n},k}P_{\vec{n}-\vec{e}_{k}}(x).$

It is straightforward to see that the compatibility condition of the Lax set
(2.3) gives us the nonlinear system containing $r^{2}+r$ equations, even
though this system apparently looks like an overdetermined system. However,
with the help of the consistency relations (1.7), the $r^{2}+r$ equations are
reduced to $2r$ equations and the evolution of this system is thus uniquely
determined. Arranging these arguments, we again arrive at the equations (1.9)
of Theorem 1.1:

$\displaystyle\frac{d}{dt}a_{\vec{n},k}=a_{\vec{n},k}[b_{\vec{n}-\vec{e}_{k},k}-b_{\vec{n},k}],$
$\displaystyle\frac{d}{dt}b_{\vec{n},k}=\sum_{j=1}^{r}(a_{\vec{n},j}-a_{\vec{n}+\vec{e}_{k},j}),\qquad
1\leq k\leq r.$ (2.18)

Recalling the determinant expression of the solution (2.13), one can easily
find that the coefficients $a_{\vec{n},j}$ and $b_{\vec{n},j}$ are also
expressed as follows

$\displaystyle
b_{\vec{n},j}=\frac{\frac{d}{dt}\tau_{\vec{n}}}{\tau_{\vec{n}}}-\frac{\frac{d}{dt}\tau_{\vec{n}+\vec{e}_{j}}}{\tau_{\vec{n}+\vec{e}_{j}}},\qquad
a_{\vec{n},j}=\frac{\tau_{\vec{n}+\vec{e}_{j}}\tau_{\vec{n}-\vec{e}_{j}}}{\tau_{\vec{n}}^{2}}.$

Substituting this expression into (2.18), we can can get the following
statement.

######  Proposition 2.2.

The $\tau$-function verifies the following bilinear equation

$\displaystyle\tau_{\vec{n}}\frac{d^{2}}{dt^{2}}\tau_{\vec{n}}=\left(\frac{d}{dt}\tau_{\vec{n}}\right)^{2}+\sum_{k=1}^{r}\tau_{\vec{n}+\vec{e}_{k}}\tau_{\vec{n}-\vec{e}_{k}}.$

This bilinear equation is also a generalization of that of the ordinary Toda
equation.

###  2.4 Toda chains for the diagonal m-OPs

Recall that $r$-OPs (or diagonal, or step-line m-OPs) could be recovered from
the lattice of m-OPs using (1.10). One can thus derive another integrable
system, especially related to $r$-OPs, from the m-Toda lattice (2.18). Here we
illustrate this taking the case $r=2$ for simplicity. Let us take in (1.5) two
different linear functionals $\mathcal{L}_{1}$, $\mathcal{L}_{2}$ and denote
the corresponding m-OPs by $P_{m,n}$. Then the recurrence relations (1.6) of
m-OPs $P_{m,n}$ take the form

$\displaystyle
xP_{m,n}(x)=P_{m+1,n}(x)+b_{m,n,1}P_{m,n}(x)+a_{m,n,1}P_{m-1,n}(x)+a_{m,n,2}P_{m,n-1}(x),$
$\displaystyle
xP_{m,n}(x)=P_{m,n+1}(x)+b_{m,n,2}P_{m,n}(x)+a_{m,n,1}P_{m-1,n}(x)+a_{m,n,2}P_{m,n-1}(x),$
(2.19)

and the corresponding integrable system (1.9) is

$\displaystyle\dot{b}_{m,n,1}=a_{m+1,n,1}-a_{m,n,1}+a_{m+1,n,2}-a_{m,n,2},$
$\displaystyle\dot{b}_{m,n,2}=a_{m,n+1,1}-a_{m,n,1}+a_{m,n+1,2}-a_{m,n,2},$
$\displaystyle\dot{a}_{m,n,1}=a_{m,n,1}(b_{m,n,2}-b_{m-1,n,2}),\qquad\dot{a}_{m,n,2}=a_{m,n,2}(b_{m,n,2}-b_{m,n-1,2}),$
(2.20)

with the contiguous relations for the initial values

$\displaystyle b_{m,n+1,1}-b_{m,n,1}=b_{m+1,n,2}-b_{m,n,2},$ $\displaystyle
a_{m+1,n,1}-a_{m,n+1,1}+a_{m+1,n,2}-a_{m,n+1,2}=b_{m,n,1}b_{m+1,n,2}-b_{m,n+1,1}b_{m,n,2},$
$\displaystyle
a_{m,n+1,1}(b_{m-1,n,1}-b_{m-1,n,2})=a_{m,n,1}(b_{m,n,1}-b_{m,n,2}),$
$\displaystyle
a_{m+1,n,1}(b_{m,n-1,1}-b_{m,n-1,2})=a_{m,n,2}(b_{m,n,1}-b_{m,n,2}).$ (2.21)

We get the 2-orthogonal polynomials $\\{q_{n}\\}$ in the following manner

$\displaystyle q_{2n}(x)=p_{n,n}(x),\qquad q_{2n+1}(x)=p_{n+1,n}(x).$

It can easily be checked that $\\{q_{n}(x)\\}$ satisfy the four-term
recurrence relation

$\displaystyle
xq_{n}(x)=q_{n+1}(x)+\alpha^{(0)}_{n}q_{n}(x)+\alpha^{(1)}_{n}q_{n-1}(x)+\alpha^{(2)}_{n}q_{n-2}(x),$
(2.22)

where the coefficients are

$\displaystyle\alpha^{(0)}_{2n}=b_{n,n,1},\qquad\alpha^{(0)}_{2n+1}=b_{n+1,n,2},$
$\displaystyle\alpha^{(1)}_{2n}=a_{n,n,1}+a_{n,n,2},\qquad\alpha^{(1)}_{2n+1}=a_{n+1,n,1}+a_{n+1,n,2},$
$\displaystyle\alpha^{(2)}_{2n}=a_{n,n,1}(b_{n-1,n-1,1}-b_{n-1,n-1,2}),$
$\displaystyle\alpha^{(2)}_{2n+1}=a_{n+1,n,2}(b_{n,n-1,1}-b_{n,n-1,2}).$
(2.23)

Taking all this into account, it is straightforward to find the spectral
transformation of the 2-orthogonal polynomials $q_{n}$ as follows

$\displaystyle\dot{q}_{n}(x)=-\alpha^{(1)}_{n}q_{n-1}(x)-\alpha^{(2)}_{n}q_{n-2}(x).$
(2.24)

The two equations (2.22) and (2.24) are exactly the Lax pair of 2-orthogonal
polynomials and then the integrable system associated with 2-orthogonal
polynomials is directly derived.

######  Theorem 2.3.

We have that the following system

$\displaystyle\frac{d}{dt}\alpha^{(0)}_{n}=\alpha^{(1)}_{n+1}-\alpha^{(1)}_{n},$
$\displaystyle\frac{d}{dt}\alpha^{(1)}_{n}=\alpha^{(1)}_{n}\big{(}\alpha^{(0)}_{n}-\alpha^{(0)}_{n-1}\big{)}+\alpha^{(2)}_{n+1}-\alpha^{(2)}_{n},$
$\displaystyle\frac{d}{dt}\alpha^{(2)}_{n}=\alpha^{(2)}_{n}\big{(}\alpha^{(0)}_{n}-\alpha^{(0)}_{n-2}\big{)},$
(2.25)

is satisfied. Moreover, one can rewrite the system in the Lax form

$\displaystyle\frac{d}{dt}L:=[L,(L)_{-}],\qquad
L:=\begin{pmatrix}\alpha^{(0)}_{0}&1&&&\\\
\alpha^{(1)}_{1}&\alpha^{(0)}_{1}&1&&\\\
\alpha^{(2)}_{2}&\alpha^{(1)}_{2}&\alpha^{(0)}_{2}&1&\\\
&\ddots&\ddots&\ddots&\ddots\end{pmatrix},$

where $(X)_{-}$ denotes the strictly lower part of the semi-infinite matrix
$X$.

This system is exactly the special case of the full Kostant–Toda lattice
investigated in [15] (see also [4]). It should be noted that the case $r=2$ is
discussed in [15], while our method is valid for the general case $r\geq 2$.

######  Remark 2.4.

The correspondence (2.23) is nothing but the Miura transformation from the
m-Toda lattice (2.20) to the special case of the full Kostant–Toda lattice
(2.25).

###  2.5 Example

We shall exhibit an interesting exact solution to the m-Toda lattice (2.20).
Let us introduce the m-OPs $P_{m,n}$ which satisfy the following multiple
orthogonality relation

$\displaystyle\mathcal{L}_{1}\big{[}x^{i}P_{m,n}(x)\big{]}=\int_{0}^{\infty}x^{i+\delta}P_{m,n}(x)e^{-x(t+\kappa_{1})}dx=0,\qquad
i=0,\dots,m-1,$
$\displaystyle\mathcal{L}_{2}\big{[}x^{j}P_{m,n}(x)\big{]}=\int_{0}^{\infty}x^{j+\delta}P_{m,n}(x)e^{-x(t+\kappa_{2})}dx=0,\qquad
j=0,\dots,n-1,$

where $\delta>-1$, $\kappa_{1}\neq\kappa_{2}$ and $t>-\kappa_{i}$ is assumed
for $i=1,2$. The corresponding m-OPs belong to the class of multiple Laguerre
polynomials of the second kind [12] and they are shown [55] to satisfy the
nearest-neighbor recurrence relation (2.19) with

$\displaystyle
b_{m,n,1}=\frac{2m+n+\delta+1}{\kappa_{1}+t}+\frac{n}{\kappa_{2}+t},\qquad
b_{m,n,2}=\frac{m+2n+\delta+1}{\kappa_{2}+t}+\frac{m}{\kappa_{1}+t},$
$\displaystyle a_{m,n,1}=\frac{m(m+n+\delta)}{(\kappa_{1}+t)^{2}},\qquad
a_{m,n,2}=\frac{n(m+n+\delta)}{(\kappa_{2}+t)^{2}}.$ (2.26)

It is easily verified that the corresponding linear functionals
$\mathcal{L}_{1}$, $\mathcal{L}_{2}$ have the property (2.15). Hence the
coefficients (2.26) directly give the special solutions to the m-Toda lattice
(2.20). Furthermore, by using (2.23), we can also obtain the corresponding
solution to the special case of the full Kostant–Toda lattice (2.25) as
follows

$\displaystyle\alpha^{(0)}_{2n}=\frac{3n+\delta+1}{\kappa_{1}+t}+\frac{n}{\kappa_{2}+t},\qquad\alpha^{(0)}_{2n+1}=\frac{3n+\delta+2}{\kappa_{2}+t}+\frac{n+1}{\kappa_{1}+t},$
$\displaystyle\alpha^{(1)}_{2n}=\frac{n(2n+\delta)}{(\kappa_{1}+t)^{2}}+\frac{n(2n+\delta)}{(\kappa_{2}+t)^{2}},\qquad\alpha^{(1)}_{2n+1}=\frac{(n+1)(2n+\delta+1)}{(\kappa_{1}+t)^{2}}+\frac{n(2n+\delta+1)}{(\kappa_{2}+t)^{2}},$
$\displaystyle\alpha^{(2)}_{2n}=(\kappa_{2}-\kappa_{1})\frac{n(2n+\delta)(2n+\delta-1)}{(\kappa_{1}+t)^{3}(\kappa_{2}+t)},\qquad\alpha^{(2)}_{2n+1}=(\kappa_{1}-\kappa_{2})\frac{n(2n+\delta)(2n+\delta+1)}{(\kappa_{1}+t)(\kappa_{2}+t)^{3}}.$

Interestingly, it is obvious that these solutions have a pole at
$t=-\kappa_{1}$ and $t=-\kappa_{2}$, which shows the presence of singularities
at certain finite times.

##  3 The discrete-time higher analogues of the Toda lattice

###  3.1 Discrete-time Toda equations and the q-d algorithm

In this section we recast two approaches presented in [46, 49, 50, 51]. The
first method, which was presented in [46] and [50], will then be generalized
in Subsection 3.2 to the case of multiple orthogonal polynomials. As for the
second one [49], we will show in this section how to modify this to be
applicable in the case of multiple orthogonal polynomials. Finally, following
the discrete integrability approach proposed in [20], we will adapt it to the
settings in question in Subsection 3.4.

To start with, suppose we are given a positive measure $d\mu$ on $(0,+\infty)$
for which all the moments exist. Clearly, the measure $x^{t}d\mu(x)$, which is
defined on $(0,+\infty)$, is also positive for all $t\in{\mathbb{Z}}_{+}$. Let
$P_{n}^{t}$ be the family of polynomials orthogonal with respect to the
measure $x^{t}d\mu(x)$ on $(0,+\infty)$. Introducing the moments $\mu_{j}$ of
the measure $d\mu$

$\displaystyle\mu_{j}=\int_{0}^{\infty}x^{j}d\mu(x),\qquad j=0,1,\dots,$

one can easily check that the monic orthogonal polynomials $P_{n}^{t}$ can be
presented in the following manner

$\displaystyle
P_{n}^{t}(x)=\frac{1}{\tau_{n}^{t}}\left|\begin{matrix}\mu_{t}&\ldots&\mu_{n+t-1}&\mu_{n+t}\\\
\vdots&\vdots&\vdots&\vdots\\\ \mu_{n+t-1}&\ldots&\mu_{2n+t-2}&\mu_{2n+t-1}\\\
1&\ldots&x^{n-1}&x^{n}\end{matrix}\right|,$

with the corresponding Hankel determinant

$\displaystyle\tau_{n+1}^{t}=\left|\begin{matrix}\mu_{t}&\ldots&\mu_{n+t}\\\
\vdots&\vdots&\vdots\\\ \mu_{n+t}&\ldots&\mu_{2n+t}\end{matrix}\right|.$

It is well known that orthogonal polynomials are related by means of
recurrence relations. To get those relations in the form that will be used
here, let us recall the Sylvester identity

$\displaystyle|A||A_{r,s;p,q}|=|A_{r;p}||A_{s;q}|-|A_{r;q}||A_{s;p}|,$

where $|A|$ stands for the determinant of the $j\times j$ matrix $A$ and
$A_{r,s;p,q}$ denotes the submatrix of $A$ formed by deleting columns number
$r$, $s$ and rows number $p$, $q$; $A_{\alpha;\beta}$ denotes the submatrix of
$A$ that is obtained from $A$ by removing the $\alpha$th column and $\beta$th
row.

Next, applying two different forms of the Sylvester identity

$\displaystyle|A||A_{1,n+1;1,n+1}|=|A_{1;1}||A_{n+1;n+1}|-|A_{1;n+1}||A_{n+1;1}|,$
$\displaystyle|A||A_{n,n+1;1,n+1}|=|A_{n;1}||A_{n+1;n+1}|-|A_{n;n+1}||A_{n+1;n+1}|$

to the determinant $\tau_{n}^{t}P_{n}^{t}$ leads to the relations

$\displaystyle P_{n+1}^{t}(x)=xP_{n}^{t+1}(x)-V_{n}^{t}P_{n}^{t},$ (3.1)
$\displaystyle P_{n+1}^{t}(x)=xP_{n}^{t+2}(x)-W_{n}^{t}P_{n}^{t+1},$ (3.2)

where

$\displaystyle
V_{n}^{t}=\frac{\tau_{n+1}^{t+1}\tau_{n}^{t}}{\tau_{n}^{t+1}\tau_{n+1}^{t}},\qquad
W_{n}^{t}=\frac{\tau_{n+1}^{t+1}\tau_{n}^{t+1}}{\tau_{n+1}^{t}\tau_{n}^{t+2}}.$

The transformation (3.1) from $P_{n}^{t}$ to $P_{n}^{t+1}$ is called the
Christoffel transformation (for instance see [21] that has a review of the
theory of such transformations). The idea of the transformation is to
construct polynomials orthogonal with respect to $xd\mu(x)$ provided that the
polynomials orthogonal with respect to $d\mu(x)$ are given. Usually, the
Christoffel transformation appears in the context of Christoffel–Darboux
kernels rather than the Sylvester identity and that is why this transformation
is called Christoffel transformation. More precisely, $P_{n}^{t+1}$ can be
represented by means of the Christoffel–Darboux kernel in the following manner

$\displaystyle
P_{n}^{t+1}(x)=\frac{1}{P_{n}^{t}(0)}\frac{P_{n+1}^{t}(x)P_{n}^{t}(0)-P_{n+1}^{t}(0)P_{n}^{t}(x)}{x},$

which is just another form of (3.1).

The reciprocal to the Christoffel transformation is called the Geronimus
transformation and it has the following form (for instance see [21] and [51])

$\displaystyle P_{n}^{t}(x)=P_{n}^{t+1}(x)+B_{n}^{t}P_{n-1}^{t+1}(x).$ (3.3)

Setting

$\displaystyle A_{n}^{t}=-\frac{P_{n+1}^{t}(0)}{P_{n}^{t}(0)}$

we see that the consistency of the Christoffel transformation

$\displaystyle P_{n}^{t+1}(x)=\frac{P_{n+1}^{t}(x)+A_{n}^{t}P_{n}^{t}(x)}{x}$
(3.4)

and the Geronimus transformation (3.3) leads to the q-d algorithm or, which is
equivalent, to the discrete-time Toda equation (1.2) [46]. Indeed, on the one
hand substituting (3.4) to (3.3) gives

$\displaystyle
xP_{n}^{t}(x)=P_{n+1}^{t}(x)+\big{(}A_{n}^{t}+B_{n}^{t}\big{)}P_{n}^{t}(x)+A_{n-1}^{t}B_{n}^{t}P_{n-1}^{t}.$

On the other hand, if we plug (3.3) into (3.4) we get

$\displaystyle
xP_{n}^{t+1}(x)=P_{n+1}^{t+1}(x)+\big{(}A_{n}^{t}+B_{n+1}^{t}\big{)}P_{n}^{t+1}(x)+A_{n}^{t}B_{n}^{t}P_{n-1}^{t+1}.$

Now, comparing the corresponding coefficients leads to (1.2):

$\displaystyle A_{n}^{t+1}+B_{n}^{t+1}=A_{n}^{t}+B_{n+1}^{t},\qquad
A_{n-1}^{t+1}B_{n}^{t+1}=A_{n}^{t}B_{n}^{t},\qquad t,n\in{\mathbb{Z}}_{+}.$

In fact, this idea is the first scheme that we are going to extend to the
setting of multiple orthogonal polynomials in the next subsection.

The second approach we mention here was presented in [49] for orthogonal
polynomials. The authors propose to use the transformation (3.2) rather than
(3.3) and they show that the relations (3.1) and (3.2) give a Lax pair for the
discretization of the Toda chain. Roughly speaking, the consistency of (3.1)
and (3.2) leads to a discrete zero curvature condition (for more information
about discrete integrability and zero curvature conditions see [3, 20, 46,
49]).

Now let us quickly see how it works. At first, introduce the wave function

$\displaystyle\widetilde{\Psi}_{n,t}(x)=\big{(}P^{t}_{n}(x),P^{t+1}_{n}(x)\big{)}^{\top}.$

Next, using (3.1) and (3.2) we derive

$\displaystyle\widetilde{\Psi}_{n+1,t}=\widetilde{L}_{n,t}\widetilde{\Psi}_{n,t},\qquad\widetilde{\Psi}_{n,t+1}=\widetilde{M}_{n,t}\widetilde{\Psi}_{n,t},$
(3.5)

where the transition matrices are defined as follows

$\displaystyle\widetilde{L}_{n,t}=\left(\begin{matrix}-V_{n}^{t}&x\\\
-V_{n}^{t}&x+W_{n}^{t}-V_{n}^{t+1}\end{matrix}\right),\qquad\widetilde{M}_{n,t}=\frac{1}{x}\left(\begin{matrix}0&x\\\
-V_{n}^{t}&x+W_{n}^{t}\end{matrix}\right).$

The consistency of the linear systems (3.5) is then equivalent to the zero
curvature condition

$\displaystyle
0=\widetilde{L}_{n,t+1}\widetilde{M}_{n,t}-\widetilde{M}_{n+1,t}\widetilde{L}_{n,t},$
(3.6)

which can be simplified to the quotient-difference scheme

$\displaystyle V_{n}^{t+2}+W_{n+1}^{t}=V_{n+1}^{t}+W_{n}^{t+1},\qquad
W_{n}^{t}V_{n+1}^{t}=V_{n}^{t+1}W_{n+1}^{t},$ (3.7)

which, as was already mentioned, can be considered as the discrete-time Toda
equation [46, 49]. While on the subject, let us mention that the quotient-
difference scheme, along with many other relations between orthogonal
polynomials, naturally occur in the context of Padé tables [28].

Now we are in the position to modify the approach we just recalled. The reason
to do that is the fact that instead of having (3.1) and (3.2) one is usually
given the three-term recurrence relation

$\displaystyle
xP_{n}^{t}(x)=P_{n+1}^{t}(x)+b_{n}^{t}P_{n}^{t}(x)+a_{n}^{t}P_{n-1}^{t}(x).$
(3.8)

At the same time, it is not so hard to check that the combination of (3.1) and
(3.2) leads to the monic version of the three-term recurrence relation (see
[49])

$\displaystyle
xP_{n}^{t}=P_{n+1}^{t}+\big{(}V_{n}^{t}+V_{n-1}^{t+1}-W_{n}^{t}\big{)}P_{n}^{t}+\big{(}V_{n-1}^{t+1}-W_{n-1}^{t}\big{)}V_{n-1}^{t}P_{n-1}^{t}.$

Hence, we also have formulas for the coefficients $a_{n}^{t}$ and $b_{n}^{t}$
in terms of Hankel determinants

$\displaystyle
a_{n}^{t}=\big{(}V_{n-1}^{t+1}-W_{n-1}^{t}\big{)}V_{n-1}^{t},\qquad
b_{n}^{t}=V_{n}^{t}+V_{n-1}^{t+1}-W_{n}^{t}.$ (3.9)

Once we have the coefficients of the three-term recurrence relation, it is in
many cases a simple task to reconstruct the coefficients of the Christoffel
transformation (we can use either the determinant formula given in (3.1) or
the formula in terms of the polynomials based on (3.4)). Thus, in order to
find the Lax pair, it is preferable to use (3.8) and (3.1).

######  Proposition 3.1.

Let us consider the following vector-valued wave function

$\displaystyle{\Psi}_{n,t}(x)=\big{(}P^{t}_{n}(x),P^{t}_{n-1}(x)\big{)}^{\top}.$

Then the corresponding transition matrices are

$\displaystyle L_{n,t}=\left(\begin{matrix}b_{n}^{t}-x&a_{n}^{t}\\\
1&0\end{matrix}\right),$ (3.10) $\displaystyle
M_{n,t}=\frac{1}{x}\left(\begin{matrix}b_{n}^{t}+V_{n}^{t}-x&a_{n}^{t}\\\
1&V_{n-1}^{t}\end{matrix}\right),$ (3.11)

and they give another Lax pair for the discrete time Toda equation (3.7).

###### Proof.

To see that the statement holds, we notice that, after some manipulations with
(3.8) and (3.1), one can get the following equalities

$\displaystyle{\Psi}_{n+1,t}=L_{n,t}{\Psi}_{n,t},\qquad{\Psi}_{n,t+1}=M_{n,t}{\Psi}_{n,t},$

where the transition matrices $L_{n,t}$ and $M_{n,t}$ are given by (3.10) and
(3.11), respectively. Next, since we have the relation

$\displaystyle\widetilde{\Psi}_{n,t}(x)=\left(\begin{matrix}1&0\\\
\frac{b_{n}^{t}+V_{n}^{t}}{x}-1&\frac{a_{n}^{t}}{x}\end{matrix}\right){\Psi}_{n,t}(x),$

it is clear that

$\displaystyle 0={L}_{n,t+1}{M}_{n,t}-{M}_{n+1,t}{L}_{n,t}$

is equivalent to (3.6) and, in turn, reduces to (3.7). ∎

###  3.2 Christoffel and Geronimus transformations for m-OPs

In this section, we will use the first method from the previous subsection to
get an integrable discretization of the m-Toda lattice (2.18). As for the
discretization of integrable systems, many techniques have been proposed and
investigated (see for the details, e.g., [29, 52]). Nonetheless, it is quite
convenient to construct the discretization by means of the discrete spectral
transformation of m-OPs as was done in [46, 50, 51].

Let us work on the discrete spectral transformations of m-OPs, a mapping from
m-OPs to another m-OPs. It is not so difficult to get to a generalization of
the Christoffel transformation.

######  Proposition 3.2.

Let $\\{P_{\vec{n}}^{t}(x)\\}$ be m-OPs at some time $t$ with respect to
$\mathcal{L}_{1}^{t},\dots,\mathcal{L}_{r}^{t}$ and define the new sequence of
polynomials $\\{P_{\vec{n},j}^{t+1}(x)\\}$ for $j=1,\dots,r$ by

$\displaystyle
P_{\vec{n},j}^{t+1}(x)=\frac{1}{x-\lambda_{t}}\big{(}P_{\vec{n}+\vec{e}_{j}}^{t}(x)+A_{\vec{n},j}^{t}P_{\vec{n}}^{t}(x)\big{)},\qquad
A_{\vec{n},j}^{t}=-\frac{P_{\vec{n}+\vec{e}_{j}}^{t}(\lambda_{t})}{P_{\vec{n}}^{t}(\lambda_{t})},$
(3.12)

where $\lambda_{t}\in\\{z\in\mathbb{R}\,|\,P_{\vec{n}}^{t}(z)\neq
0,\,\textrm{for all}\,\,\vec{n}\in\mathbb{Z}_{+}^{r}\\}$. Then

$\displaystyle
P_{\vec{n},1}^{t+1}(x)=\cdots=P_{\vec{n},r}^{t+1}(x)=P_{\vec{n}}^{t+1}(x)$
(3.13)

holds and $\\{P_{\vec{n}}^{t+1}(x)\\}$ are again m-OPs with respect to the new
linear functional $\mathcal{L}_{1}^{t+1},\dots,\mathcal{L}_{r}^{t+1}$ defined
by

$\displaystyle\mathcal{L}_{j}^{t+1}[\,\cdot\,]:=\mathcal{L}_{j}^{t}[(x-\lambda_{t})\,\cdot\,],\qquad
j=1,\dots,r.$ (3.14)

###### Proof.

From the multiple orthogonality, it is easy to verify

$\displaystyle\mathcal{L}_{i}^{t+1}\big{[}x^{k}P_{\vec{n},j}^{t+1}(x)\big{]}=\mathcal{L}_{i}^{t}\big{[}x^{k}(x-\lambda_{t})P_{\vec{n},j}^{t+1}(x)\big{]}=\mathcal{L}_{i}^{t}\big{[}x^{k}\big{(}P_{\vec{n}+\vec{e}_{j}}^{t}(x)+A_{\vec{n},j}^{t}P_{\vec{n},j}^{t}(x)\big{)}\big{]}$
$\displaystyle\hphantom{\mathcal{L}_{i}^{t+1}\big{[}x^{k}P_{\vec{n},j}^{t+1}(x)\big{]}}{}=0,\qquad
k=0,1,\dots,n_{i}-1,$

for all $i,j=1,\dots,r$. Then the uniqueness of monic m-OPs gives us (3.13)
and this completes the proof. ∎

In case $r=1$, the transformation (3.12) coincides with the Christoffel
transformation for OPs. We shall refer to the transformation (3.12) as the
Christoffel transformation for m-OPs. Setting the initial time $t=0$ and
iterating the Christoffel transformation

$\displaystyle
P_{\vec{n}}^{t+1}(x)=\frac{1}{x-\lambda_{t}}\big{(}P_{\vec{n}+\vec{e}_{j}}^{t}(x)+A_{\vec{n},j}^{t}P_{\vec{n}}^{t}(x)\big{)},$
$\displaystyle
A_{\vec{n},j}^{t}=-\frac{P_{\vec{n}+\vec{e}_{j}}^{t}(\lambda_{t})}{P_{\vec{n}}^{t}(\lambda_{t})},\qquad
j=1,\dots,r,$ (3.15)

we can obtain the chain of m-OPs

$\displaystyle\big{\\{}P_{\vec{n}}^{0}(x)\big{\\}}\rightarrow\big{\\{}P_{\vec{n}}^{1}(x)\big{\\}}\rightarrow\cdots\rightarrow\big{\\{}P_{\vec{n}}^{t}(x)\big{\\}}\rightarrow\big{\\{}P_{\vec{n}}^{t+1}(x)\big{\\}}\rightarrow\cdots.$

######  Remark 3.3.

We can derive the nonlinear equations from the compatibility condition of the
relations (3.15) themselves

$\displaystyle
A_{\vec{n},i}^{t+1}-A_{\vec{n},j}^{t+1}=A_{\vec{n}+\vec{e_{j}},i}^{t}-A_{\vec{n}+\vec{e}_{i},j}^{t},\qquad
A_{\vec{n},i}^{t}A_{\vec{n}+\vec{e}_{i},j}^{t}=A_{\vec{n},j}^{t}A_{\vec{n}+\vec{e}_{j},i}^{t}.$
(3.16)

As for the discrete time m-Toda lattice, the Geronimus transformation, which
in a way is the reciprocal to the Christoffel transformation, plays a key
role. The Geronimus transformation for m-OPs is given in the following
statement.

######  Theorem 3.4.

The sequence of m-OPs $\\{P_{\vec{n}}^{t}(x)\\}$ connected by the relation
(3.15) satisfies the following contiguous relations

$\displaystyle
P_{\vec{n}}^{t}(x)=P_{\vec{n}}^{t+1}(x)+\sum_{k=1}^{r}B_{\vec{n},k}^{t}P_{\vec{n}-\vec{e}_{k}}^{t+1}(x),\qquad
B_{\vec{n},k}^{t}=-a_{\vec{n},k}^{t}\frac{P_{\vec{n}-\vec{e}_{k}}^{t}(\lambda_{t})}{P_{\vec{n}}^{t}(\lambda_{t})},$
(3.17)

where $a_{\vec{n},k}^{t}$ are the coefficients of the nearest neighbor
recurrence relation

$\displaystyle
xP_{\vec{n}}^{t}(x)=P_{\vec{n}+\vec{e}_{j}}^{t}(x)+b_{\vec{n},j}^{t}P_{\vec{n}}^{t}(x)+\sum_{k=1}^{r}a_{\vec{n},k}^{t}P_{\vec{n}-\vec{e}_{k}}^{t}(x),\qquad
j=1,\dots,r.$ (3.18)

###### Proof.

First, consider the polynomial $P_{\vec{n}}^{t}(x)-P_{\vec{n}}^{t+1}(x)$,
which is a polynomial of degree $|\vec{n}|-1$. From the multiple orthogonality
relation and (3.14), we can easily check

$\displaystyle\mathcal{L}_{i}^{t+1}\big{[}x^{j}P_{\vec{n}}^{t}(x)\big{]}=\mathcal{L}_{i}^{t}\big{[}(x-\lambda_{t})x^{j}P_{\vec{n}}^{t}(x)\big{]}=0,\qquad
j=0,\dots,n_{i}-2.$ (3.19)

This readily shows $P_{\vec{n}}^{t}(x)-P_{\vec{n}}^{t+1}(x)$ is orthogonal to
all polynomials of degree less than $n_{j}-1$ with respect to the linear
functional $\mathcal{L}_{j}^{t+1}$. Hence, we can write
$P_{\vec{n}}^{t}(x)-P_{\vec{n}}^{t+1}(x)$ as a linear combination of the
polynomials $P^{t+1}_{\vec{n}-\vec{e}_{j}}$, $j=1,\dots,r$, which form a basis
for the linear space of all polynomials of degree less than $|\vec{n}|$ and
satisfy the multiple orthogonality conditions

$\displaystyle\mathcal{L}_{j}^{t+1}\big{[}x^{k}P_{\vec{n}}^{t+1}(x)\big{]}=0,\qquad
j=1,\dots,r,$

for $k\leq n_{j}-2$. We write
$P_{\vec{n}}^{t}(x)-P_{\vec{n}}^{t+1}(x)=\sum\limits_{k=1}^{r}B_{\vec{n},k}^{t}P_{\vec{n}-\vec{e}_{k}}^{t+1}(x)$.
From (3.19), we can easily find

$\displaystyle
B_{\vec{n},k}^{t}=\frac{\mathcal{L}_{k}^{t+1}[x^{n_{k}-1}P_{\vec{n}}^{t}(x)]}{\mathcal{L}_{k}^{t+1}[x^{n_{k}-1}P_{\vec{n}-\vec{e}_{k}}^{t+1}(x)]},\qquad
k=1,\dots,r.$ (3.20)

Here, some calculations show

$\displaystyle\mathcal{L}_{k}^{t+1}\big{[}x^{n_{k}-1}P_{\vec{n}}^{t}(x)\big{]}=\mathcal{L}_{k}^{t}\big{[}(x-\lambda_{t})x^{n_{k}-1}P_{\vec{n}}^{t}(x)\big{]}=\mathcal{L}_{k}^{t}\big{[}x^{n_{k}}P_{\vec{n}}^{t}(x)\big{]},$
(3.21)

and

$\displaystyle\mathcal{L}_{k}^{t+1}\big{[}x^{n_{k}-1}P_{\vec{n}-\vec{e}_{k}}^{t+1}(x)\big{]}=\mathcal{L}_{k}^{t}\big{[}(x-\lambda_{t})x^{n_{k}-1}P_{\vec{n}-\vec{e}_{k}}^{t+1}(x)\big{]}=\mathcal{L}_{k}^{t}\big{[}x^{n_{k}-1}\big{(}P_{\vec{n}}^{t}(x)+A_{\vec{n},k}^{t}P_{\vec{n}-\vec{e}_{k}}^{t}(x)\big{)}\big{]}$
$\displaystyle\hphantom{\mathcal{L}_{k}^{t+1}\big{[}x^{n_{k}-1}P_{\vec{n}-\vec{e}_{k}}^{t+1}(x)\big{]}}{}=A_{\vec{n},k}^{t}\mathcal{L}_{k}^{t}\big{[}x^{n_{k}-1}P_{\vec{n}-\vec{e}_{k}}^{t}(x)\big{]}.$
(3.22)

We can also calculate the coefficients of (3.18) from the multiple
orthogonality relation

$\displaystyle
a_{\vec{n},k}^{t}=\frac{\mathcal{L}_{k}^{t}[x^{n_{k}}P_{\vec{n}}^{t}(x)]}{\mathcal{L}_{k}^{t}[x^{n_{k}-1}P_{\vec{n}-\vec{e}_{k}}^{t}(x)]}.$
(3.23)

Combining (3.20)–(3.23), we finally arrive at the following result:

$\displaystyle
B_{\vec{n}}^{t}=\frac{a_{\vec{n},k}^{t}}{A_{\vec{n},k}^{t}}=-a_{\vec{n},k}^{t}\frac{P_{\vec{n}-\vec{e}_{k}}^{t}(\lambda_{t})}{P_{\vec{n}}^{t}(\lambda_{t})},\qquad
k=1,\dots,r.$

This completes the proof. ∎

######  Remark 3.5.

Although we can formally derive all the relations in Section 3 for
$\lambda_{t}$ that changes with the discrete time $t$, we are only concerned
with the case $\lambda_{t}=\lambda$ and, particularly, $\lambda_{t}=0$. The
reason is that the basis for our construction of integrable systems is a
perfect system of measures $d\mu_{1},\dots,d\mu_{r}$ and it has to be perfect
at all times. In other words, we need to have all the multiple orthogonal
polynomials to exist for any multi-index for all values of the time. However,
it is still an open question when the system remains perfect under Christoffel
or Geronimus transformations. Moreover, the only examples we know at the
moment correspond to the case $\lambda_{t}=0$, which can be easily modified to
$\lambda_{t}=\lambda$ (see Subsection 3.5).

From (3.15) and (3.17), we can reproduce the nearest neighbor recurrence
relation (3.18) and the coefficients can explicitly be written in terms of
$\big{\\{}A_{\vec{n},j}^{t},B_{\vec{n},j}^{t}\big{\\}}$

$\displaystyle
b_{\vec{n},j}^{t}=A_{\vec{n},j}^{t}+\sum_{k=1}^{r}B_{\vec{n},k}^{t}+\lambda_{t},\qquad
a_{\vec{n},j}^{t}=A_{\vec{n}-\vec{e}_{j},j}^{t}B_{\vec{n},j}^{t},\qquad
j=1,\dots,r.$ (3.24)

Substituting (3.24) into (1.7) and also using the relation (3.16), we obtain,
after some calculations and simplifications, the contiguous relations for
$\\{A_{\vec{n},j}^{t},B_{\vec{n},j}^{t}\\}$.

######  Corollary 3.6.

If we put $\lambda_{t}=0$, the coefficients
$\big{\\{}A_{\vec{n},j}^{t},B_{\vec{n},j}^{t}\big{\\}}$ in (3.15) and (3.17)
satisfy the following difference equations on $\vec{n}$

$\displaystyle
A_{\vec{n},i}^{t}A_{\vec{n}+\vec{e}_{i},j}^{t}=A_{\vec{n},j}^{t}A_{\vec{n}+\vec{e}_{j},i}^{t},$
$\displaystyle
A_{\vec{n}+\vec{e}_{j},i}^{t}-A_{\vec{n}+\vec{e_{i}},j}^{t}+A_{\vec{n},j}^{t}-A_{\vec{n},i}^{t}=\sum_{k=1}^{r}\big{(}B_{\vec{n}+\vec{e}_{i},k}^{t}-B_{\vec{n}+\vec{e}_{j},k}^{t}\big{)},$
$\displaystyle\frac{B_{\vec{n},i}^{t}}{B_{\vec{n}+\vec{e}_{j},i}^{t}}=\frac{A_{\vec{n},j}^{t}-A_{\vec{n}+\vec{e}_{j}-\vec{e}_{i},i}^{t}}{A_{\vec{n},j}^{t}-A_{\vec{n},i}^{t}}$
(3.25)

for all $t$ and $i,j=1,\dots,r$.

From the Christoffel and Geronimus transformations for m-OPs, we obtain the
discrete Lax set for which the m-OPs appear as their eigenfunctions

$\displaystyle(x-\lambda_{t})P_{\vec{n}}^{t+1}(x)=P_{\vec{n}+\vec{e}_{j}}^{t}(x)+A_{\vec{n},j}^{t}P_{\vec{n}}^{t}(x),\qquad
j=1,\dots,r,$ $\displaystyle
P_{\vec{n}}^{t}(x)=P_{\vec{n}}^{t+1}(x)+\sum_{k=1}^{r}B_{\vec{n},k}^{t}P_{\vec{n}-\vec{e}_{k}}^{t+1}(x).$
(3.26)

Next, the compatibility condition for the relations in (3.26) gives us
$r^{2}+r$ equations and, therefore, the obtained system seems overdetermined.
However, if we take into account the contiguous relation (3.25), then these
$r^{2}+r$ equations are reduced to $2r$ equations and the evolution is
uniquely determined. Summing up these arguments, we get to Theorem 1.5:

$\displaystyle
A_{\vec{n},j}^{t+1}+\sum_{k=1}^{r}B_{\vec{n},k}^{t+1}=A_{\vec{n},j}^{t}+\sum_{k=1}^{r}B_{\vec{n}+\vec{e_{j}},k}^{t},$
$\displaystyle
A_{\vec{n}-\vec{e}_{j},j}^{t+1}B_{\vec{n},j}^{t+1}=A_{\vec{n},j}^{t}B_{\vec{n},j}^{t},\qquad
1\leq j\leq r.$ (3.27)

######  Remark 3.7.

As with the continuous time m-Toda lattice, we see here that the system (3.27)
is solvable (integrable) only if the initial values are chosen so that (3.25)
are satisfied.

One can also see that in the case $r=1$ the system (3.27) coincides with the
discrete time Toda lattice (1.2). Therefore, it is natural to call (3.27) a
discrete multiple Toda (dm-Toda) lattice.

######  Remark 3.8.

The dm-Toda lattice (3.27) is exactly the discrete analogue of the continuous
time m-Toda lattice (2.18). Indeed, let us introduce the new variables
$a_{\vec{n},j}(t)$ and $b_{\vec{n},j}(t)$ by the following relation

$\displaystyle\lambda_{t}=\frac{1}{\delta},\qquad
A_{\vec{n},j}^{t}=b_{\vec{n},j}(t\delta)-\frac{1}{\delta},\qquad
B_{\vec{n},j}^{t}=\delta a_{\vec{n},j}(t\delta).$ (3.28)

Then, substituting (3.28) into (3.27), we can obtain the following equations

$\displaystyle\frac{b_{\vec{n},j}(t\delta+\delta)-b_{\vec{n},j}(t\delta)}{\delta}=\sum_{k=1}^{r}a_{\vec{n}+\vec{e}_{j},k}(t\delta)-a_{\vec{n},k}(t\delta+\delta),$
$\displaystyle\frac{a_{\vec{n},j}(t\delta+\delta)-a_{\vec{n},j}(t\delta)}{\delta}=b_{\vec{n},j}(t\delta)a_{\vec{n},j}(t\delta)-b_{\vec{n}-\vec{e}_{j},j}(t\delta+\delta)a_{\vec{n},j}(t\delta+\delta).$
(3.29)

If we take $t\delta\rightarrow t$ and the continuous limit $\delta\rightarrow
0$, it is straightforward to see that the equations (3.29) go to the m-Toda
lattice (2.18).

######  Remark 3.9.

We can also verify that the contiguous relation (3.25) reduce to (1.7) in the
continuous limit after some careful calculations and simplifications.

In the previous section, we have seen that the m-Toda lattice (2.18) admits
the determinant solution (2.13) with the dispersion relation (2.14). We shall
now give the determinant solution to the dm-Toda lattice (3.27). Let us
introduce the $\tau$-function $\tau_{\vec{n}}^{t}$ defined by

$\displaystyle\tau_{\vec{n}}^{t}:=\begin{vmatrix}\mu_{0,1}^{t}&\mu_{1,1}^{t}&\cdots&\mu_{n_{1}-1,1}^{t}&\cdots&\mu_{0,r}^{t}&\cdots&\mu_{n_{r}-1,r}^{t}\\\
\mu_{1,1}^{t}&\mu_{2,1}^{t}&\cdots&\mu_{n_{1}+1,1}^{t}&\cdots&\mu_{1,r}^{t}&\cdots&\mu_{n_{r},r}^{t}\\\
\vdots&\vdots&\cdots&\vdots&\cdots&\vdots&\cdots&\vdots\\\
\mu_{|\vec{n}|-1,1}^{t}&\mu_{|\vec{n}|,1}^{t}&\cdots&\mu_{|\vec{n}|+n_{1}-2,1}^{t}&\cdots&\mu_{|\vec{n}|-1,r}^{t}&\cdots&\mu_{|\vec{n}|+n_{r}-2,r}^{t}\end{vmatrix},$

where $\mu_{i,j}^{t}:=\mathcal{L}_{j}^{t}[x^{i}]$. From (3.14), we get the
following relation

$\displaystyle\mu_{i,j}^{t+1}=\mu_{i+1,j}^{t}-\lambda_{t}\mu_{i,j}^{t},\qquad
j=1,\dots,r.$

From the determinant expression of the m-OPs (2.12) we can get by means of
elementary transformations of determinants

$\displaystyle
P_{\vec{n}}^{t}(\lambda_{t})=\frac{1}{\tau_{\vec{n}}^{t}}\begin{vmatrix}\mu_{0,1}^{t}&\cdots&\mu_{n_{1}-1,1}^{t}&\cdots&\mu_{0,r}^{t}&\cdots&\mu_{n_{r}-1,r}^{t}&1\\\
\mu_{0,1}^{t+1}&\cdots&\mu_{n_{1}-1,1}^{t+1}&\cdots&\mu_{0,r}^{t}&\cdots&\mu_{n_{r}-1,r}^{t+1}&0\\\
\mu_{1,1}^{t+1}&\cdots&\mu_{n_{1},1}^{t+1}&\cdots&\mu_{1,r}^{t}&\cdots&\mu_{n_{r},r}^{t+1}&0\\\
\vdots&\cdots&\vdots&\cdots&\vdots&\cdots&\vdots&\vdots\\\
\mu_{|\vec{n}|-1,1}^{t+1}&\cdots&\mu_{|\vec{n}|+n_{1}-2,1}^{t+1}&\cdots&\mu_{|\vec{n}|-1,r}^{t+1}&\cdots&\mu_{|\vec{n}|+n_{r}-2,r}^{t+1}&0\end{vmatrix}$
$\displaystyle\hphantom{P_{\vec{n}}^{t}(\lambda_{t})}{}=(-1)^{|\vec{n}|+1}\frac{\tau_{\vec{n}}^{t+1}}{\tau_{\vec{n}}^{t}}.$
(3.30)

Combining (2.13), (3.12), (3.17) and (3.30), we arrive at the determinant
expression of the solution:

$\displaystyle
A_{\vec{n},j}^{t}=\frac{\tau_{\vec{n}+\vec{e}_{j}}^{t+1}\tau_{\vec{n}}^{t}}{\tau_{\vec{n}+\vec{e}_{j}}^{t}\tau_{\vec{n}}^{t+1}},\qquad
B_{\vec{n},j}^{t}=\frac{\tau_{\vec{n}-\vec{e}_{j}}^{t+1}\tau_{\vec{n}+\vec{e}_{j}}^{t}}{\tau_{\vec{n}}^{t}\tau_{\vec{n}}^{t+1}},\qquad
j=1,\dots,r.$ (3.31)

We shall consider the bilinear equations for dm-Toda lattice (3.27) for the
case $\lambda_{t}=\lambda$. From (3.16), the dependent variable
$A_{\vec{n},j}^{t}$ obeys a discrete KP equation, which reduces to the
following Hirota–Miwa equation

$\displaystyle\tau_{\vec{n}+\vec{e}_{i}+\vec{e}_{j}}^{t}\tau_{\vec{n}}^{t+1}-\tau_{\vec{n}+\vec{e}_{i}}^{t}\tau_{\vec{n}+\vec{e}_{j}}^{t+1}+\tau_{\vec{n}+\vec{e}_{j}}^{t}\tau_{\vec{n}+\vec{e}_{i}}^{t+1}=0,\qquad
i\neq j.$ (3.32)

Using the Hirota–Miwa equation and substituting (3.31) into (3.27), we obtain
another bilinear equation of $\tau_{\vec{n}}^{t}$:

$\displaystyle\tau_{\vec{n}}^{t+1}\tau_{\vec{n}}^{t-1}=\big{(}\tau_{\vec{n}}^{t}\big{)}^{2}+\sum_{k=1}^{r}\tau_{\vec{n}+\vec{e}_{k}}^{t-1}\tau_{\vec{n}-\vec{e}_{k}}^{t+1},$
(3.33)

which is the multiple generalization of the bilinear equation of the ordinary
discrete Toda lattice. Indeed we can derive the dm-Toda lattice (3.27) from
these bilinear equations. Using (3.33), one has

$\displaystyle\left(\frac{\tau_{\vec{n}+\vec{e}_{j}}^{t+1}}{\tau_{\vec{n}}^{t+1}}\right)^{2}=\frac{\tau_{\vec{n}+\vec{e}_{j}}^{t+2}\tau_{\vec{n}+\vec{e}_{j}}^{t}-\tau_{\vec{n}+2\vec{e}_{j}}^{t}\tau_{\vec{n}}^{t+2}-\sum\limits_{k\neq
j}\tau_{\vec{n}+\vec{e}_{j}+\vec{e}_{k}}^{t}\tau_{\vec{n}+\vec{e}_{j}-\vec{e}_{k}}^{t+2}}{\tau_{\vec{n}}^{t+2}\tau_{\vec{n}}^{t}-\tau_{\vec{n}+\vec{e}_{j}}^{t}\tau_{\vec{n}-\vec{e}_{j}}^{t+2}-\sum\limits_{k\neq
j}\tau_{\vec{n}+\vec{e}_{k}}^{t}\tau_{\vec{n}-\vec{e}_{k}}^{t+2}},$

which can be rewritten as

$\displaystyle A_{\vec{n},j}^{t+1}+B_{\vec{n},j}^{t+1}+\sum_{k\neq
j}\frac{\tau_{\vec{n}+\vec{e_{k}}}^{t}\tau_{\vec{n}+\vec{e}_{j}}^{t+1}}{\tau_{\vec{n}+\vec{e}_{j}}^{t}\tau_{\vec{n}+\vec{e}_{k}}^{t+1}}B_{\vec{n},k}^{t+1}=A_{\vec{n},j}^{t}+B_{\vec{n}+\vec{e}_{j},j}^{t}+\sum_{k\neq
j}\frac{\tau_{\vec{n}}^{t+1}\tau_{\vec{n}+\vec{e}_{j}-\vec{e}_{k}}^{t+2}}{\tau_{\vec{n}}^{t+2}\tau_{\vec{n}+\vec{e}_{j}-\vec{e}_{k}}^{t+1}}B_{\vec{n}+\vec{e}_{j},k}^{t}.$

For all $k\neq j$, one also has

$\displaystyle
B_{\vec{n},k}^{t+1}-\frac{\tau_{\vec{n}+\vec{e_{k}}}^{t}\tau_{\vec{n}+\vec{e}_{j}}^{t+1}}{\tau_{\vec{n}+\vec{e}_{j}}^{t}\tau_{\vec{n}+\vec{e}_{k}}^{t+1}}B_{\vec{n},k}^{t+1}=B_{\vec{n}+\vec{e}_{j},k}^{t}-\frac{\tau_{\vec{n}}^{t+1}\tau_{\vec{n}+\vec{e}_{j}-\vec{e}_{k}}^{t+2}}{\tau_{\vec{n}}^{t+2}\tau_{\vec{n}+\vec{e}_{j}-\vec{e}_{k}}^{t+1}}B_{\vec{n}+\vec{e}_{j},k}^{t},$

which is easily verified from (3.32) and (3.31). Comparing the above two
relations we thus have the dm-Toda lattice (3.27). In other words, the dm-Toda
lattice is a consequence of the two bilinear equations (3.32) and (3.33).

###  3.3 The discrete-time integrable system and diagonal m-OP

Here we will introduce the Miura transformation from the dm-Toda lattice to
the discrete integrable system associated with $r$-orthogonal polynomials, as
was done in the previous section for the continuous-time m-Toda lattice. For
simplicity, we only consider the case $r=2$ and the autonomous case
$\lambda_{t}=\lambda$. We denote the corresponding m-OPs by
$\\{P_{m,n}^{t}\\}$ and we write the discrete Lax set of m-OPs as follows

$\displaystyle(x-\lambda_{t})P_{m,n}^{t+1}(x)=P_{m+1,n}^{t}(x)+A_{m,n,1}^{t}P_{m,n}^{t}(x),$
$\displaystyle(x-\lambda_{t})P_{m,n}^{t+1}(x)=P_{m,n+1}^{t}(x)+A_{m,n,2}^{t}P_{m,n}^{t}(x),$
$\displaystyle
P_{m,n}^{t}(x)=P_{m,n}^{t+1}(x)+B_{m,n,1}^{t}P_{m-1,n}^{t+1}(x)+B_{m,n,2}^{t}P_{m,n-1}^{t+1}(x).$
(3.34)

Thus the corresponding dm-Toda equation can be represented by

$\displaystyle
A_{m,n,1}^{t}+B_{m+1,n,1}^{t}+B_{m+1,n,2}^{t}=A_{m,n,1}^{t+1}+B_{m,n,1}^{t+1}+B_{m,n,2}^{t+1},$
$\displaystyle
A_{m,n,2}^{t}+B_{m,n+1,1}^{t}+B_{m,n+1,2}^{t}=A_{m,n,2}^{t+1}+B_{m,n,1}^{t+1}+B_{m,n,2}^{t+1},$
$\displaystyle
A_{m,n,1}^{t}B_{m,n,1}^{t}=A_{m-1,n,1}^{t+1}B_{m,n,1}^{t+1},\qquad
A_{m,n,2}^{t}B_{m,n,2}^{t}=A_{m,n-1,2}^{t+1}B_{m,n,2}^{t+1},$ (3.35)

with the contiguous relations of the initial values

$\displaystyle A_{m,n,1}^{0}A_{m+1,n,2}^{0}=A_{m,n+1,1}^{0}A_{m,n,2}^{0},$
$\displaystyle
A_{m,n+1,1}^{0}-A_{m,n,1}^{0}-A_{m+1,n,2}^{0}+A_{m,n,2}^{0}=B_{m+1,n,1}^{0}-B_{m,n+1,1}^{0}+B_{m+1,n,2}^{0}-B_{m,n+1,2}^{0},$
$\displaystyle\frac{B_{m,n,1}^{0}}{B_{m,n+1,1}^{0}}=\frac{A_{m-1,n+1,1}^{0}-A_{m,n,2}^{0}}{A_{m,n,1}^{0}-A_{m,n,2}^{0}},\qquad\frac{B_{m,n,2}^{0}}{B_{m+1,n,2}^{0}}=\frac{A_{m,n,1}^{0}-A_{m+1,n-1,2}^{0}}{A_{m,n,1}^{0}-A_{m,n,2}^{0}}.$

Let us introduce the sequence of $2$-orthogonal polynomials $\\{q_{n}^{t}\\}$
by the correspondence

$\displaystyle q_{2n}^{t}(x)=P_{n,n}^{t}(x),\qquad
q_{2n+1}^{t}(x)=P_{n+1,n}^{t}(x).$

In a fashion similar to the previous section, we can also directly obtain the
discrete spectral transformation of 2-orthogonal polynomials as follows

$\displaystyle(x-\lambda_{t})q_{n}^{t+1}(x)=q_{n+1}^{t}(x)+X_{n}^{t}q_{m,n}^{t}(x),\qquad
q_{n}^{t}(x)=q_{n}^{t+1}(x)+Y_{n}^{t}q_{n-1}^{t+1}(x)+Z_{n}^{t}q_{n-2}^{t+1}(x),$

with

$\displaystyle X_{2n}^{t}=A_{n,n,1}^{t},\qquad X_{2n+1}^{t}=A_{n+1,n,2}^{t},$
$\displaystyle Y_{2n}^{t}=B_{n,n,1}^{t}+B_{n,n,2}^{t},\qquad
Y_{2n+1}^{t}=B_{n+1,n,1}^{t}+B_{n+1,n,2}^{t},$ $\displaystyle
Z_{2n}^{t}=B_{n,n,1}^{t}\big{(}A_{n-1,n,1}^{t}-A_{n,n-1,2}^{t}\big{)},\qquad
Z_{2n+1}^{t}=B_{n+1,n,2}^{t}\big{(}B_{n+1,n-1}^{t}-A_{n,n,1}^{t}\big{)}.$
(3.36)

From the compatibility condition of (3.34) it follows that the discrete
integrable system associated with 2-orthogonal polynomials is given by

$\displaystyle
X_{n}^{t}+Y_{n+1}^{t}+\lambda_{t}=X_{n}^{t+1}+Y_{n}^{t+1}+\lambda_{t+1},$
$\displaystyle
X_{n}^{t}Y_{n}^{t}+Z_{n+1}^{t}=X_{n-1}^{t+1}Y_{n}^{t+1}+Z_{n}^{t+1},\qquad
X_{n}^{t}Z_{n}^{t}=X_{n-2}^{t+1}Z_{n}^{t+1}.$ (3.37)

It is easy to verify that the equation (3.37) is the integrable discretization
of the special case of Kostant–Toda equation (2.25) and the Miura
transformation from (3.35) to (3.37) is explicitly given by (3.36).

###  3.4 The consistency approach: the stationary equations
and the discrete-time dynamics

To simplify formulas and statements we restrict ourselves here to the case
where the multiple orthogonal polynomials are generated by two measures.
Nevertheless, it can straightforwardly be generalized to the case $r>2$.

Once again, recall that multiple orthogonal polynomials are a generalization
of orthogonal polynomials where the polynomials are required to be
simultaneously orthogonal with respect to two given measures [32, 55]. Now we
consider a multi-index $(n,m)\in{\mathbb{Z}}^{2}_{+}$ and suppose that
$\mu_{1}$, $\mu_{2}$ are given positive measures on the real line. Then, the
type II multiple orthogonal polynomial is the monic polynomial
$P_{n,m}(x)=x^{n+m}+\cdots$ of degree $n+m$ for which

$\displaystyle\int P_{n,m}(x)x^{j}d\mu_{1}(x)=0,\qquad j=0,1,\ldots,n-1,$
$\displaystyle\int P_{n,m}(x)x^{j}d\mu_{2}(x)=0,\qquad j=0,1,\ldots,m-1.$

As in the case of ordinary orthogonal polynomials, one can introduce the
moments

$\displaystyle\mu_{j,i}=\int x^{j}d\mu_{i}(x),\qquad i=1,2,$

and the determinant of the moment matrix

$\displaystyle\tau_{n,m}=\left|\begin{matrix}\mu_{0,1}&\cdots&\mu_{n-1,1}\\\
\mu_{1,1}&\cdots&\mu_{n,1}\\\ \vdots&\cdots&\vdots\\\
\mu_{n+m-1,1}&\cdots&\mu_{2n+m-2,1}\end{matrix}\
\begin{matrix}\mu_{0,2}&\cdots&\mu_{m-1,2}\\\ \mu_{1,2}&\cdots&\mu_{m,2}\\\
\vdots&\cdots&\vdots\\\
\mu_{n+m-1,2}&\cdots&\mu_{n+2m-2,2}\end{matrix}\right|.$

Now we see that the type II multiple orthogonal polynomial can be written as

$\displaystyle
P_{n,m}(x)=\frac{1}{\tau_{n,m}}\left|\begin{matrix}\mu_{0,1}&\cdots&\mu_{n-1,1}\\\
\mu_{1,1}&\cdots&\mu_{n,1}\\\ \vdots&\vdots&\vdots\\\
\mu_{n+m,1}&\cdots&\mu_{2n+m-1,1}\end{matrix}\
\begin{matrix}\mu_{0,2}&\cdots&\mu_{m-1,2}\\\ \mu_{1,2}&\cdots&\mu_{m,2}\\\
\vdots&\cdots&\vdots\\\ \mu_{n+m,2}&\cdots&\mu_{n+2m-1,2}\end{matrix}\
\begin{matrix}1\\\ x\\\ \vdots\\\ x^{n+m}\end{matrix}\right|$

provided that $\tau_{n,m}$ is nonvanishing. In the latter case the index
$(n,m)$ is normal. We assume that all multi-indices are normal.

In the case of multiple orthogonal polynomials, the three-term recurrence
relations are replaced with the following relation for the nearest neighbors

$\displaystyle
P_{n+1,m}(x)=(x-b_{n,m,1})P_{n,m}(x)-a_{n,m,1}P_{n-1,m}(x)-a_{n,m,2}P_{n,m-1}(x),$
$\displaystyle
P_{n,m+1}(x)=(x-b_{n,m,2})P_{n,m}(x)-a_{n,m,1}P_{n-1,m}(x)-a_{n,m,2}P_{n,m-1}(x),$
(3.38)

with $a_{0,m,1}=0$ and $a_{n,0,2}=0$ for all $n,m\geq 0$.

Unlike the case of ordinary orthogonal polynomials, the coefficients of the
recurrence relations (3.38) are solutions of a discrete integrable system even
without introducing the discrete time evolution. Here we follow the concept of
discrete integrability given in [20].

######  Proposition 3.10 ([13, 14]).

Let us consider the following vector-valued wave function

$\displaystyle\Psi_{n,m}(x)=\big{(}P_{n,m}(x),P_{n-1,m}(x),P_{n,m-1}(x)\big{)}^{\top}.$

Then the corresponding transition matrices are

$\displaystyle L_{n,m}=\begin{pmatrix}z-b_{n,m,1}&-a_{n,m,1}&-a_{n,m,2}\\\
1&0&0\\\ 1&0&b_{n,m-1,2}-b_{n,m-1,1}\end{pmatrix}$

and

$\displaystyle M_{n,m}=\begin{pmatrix}z-b_{n,m,2}&-a_{n,m,1}&-a_{n,m,2}\\\
1&b_{n-1,m,1}-b_{n-1,m,2}&0\\\ 1&0&0\end{pmatrix},$

and they give the non-trivial zero curvature condition

$\displaystyle 0=L_{n,m+1}M_{n,m}-M_{n+1,m}L_{n,m}.$ (3.39)

Indeed, it follows from (3.38) that

$\displaystyle\Psi_{n+1,m}=L_{n,m}\Psi_{n,m},\qquad\Psi_{n,m+1}=M_{n,m}\Psi_{n,m}.$
(3.40)

It is now clear that the consistency of (3.40) gives (3.39), which is in fact
a discrete integrable system [13, 14]. Namely, in [14] and [55] it is shown
that the discrete zero curvature condition (3.39) is equivalent to the
nonlinear system of difference equations (1.7) for the coefficients of the
recurrence relations (3.38). Furthermore we have the following formulas for
the recurrence coefficients [55]

$\displaystyle
a_{n,m,1}=\frac{\tau_{n+1,m}\tau_{n-1,m}}{\big{(}\tau_{n,m}\big{)}^{2}},\qquad
a_{n,m,2}=\frac{\tau_{n,m+1}\tau_{n,m-1}}{\big{(}\tau_{n,m}\big{)}^{2}},$
$\displaystyle
b_{n,m,2}-b_{n,m,1}=\frac{\tau_{n,m}\tau_{n+1,m+1}}{\tau_{n+1,m}\tau_{n,m+1}}.$
(3.41)

Nevertheless, they do not determine the coefficients of the recurrence
relations (3.38) from the moments of the given measures. Still, this obstacle
can easily be overcome.

######  Proposition 3.11.

We have that

$\displaystyle
b_{n,m+1,1}=b_{n,0,1}+\sum_{i=1}^{m}\frac{(a_{n+1,i,1}+a_{n+1,i,2})-(a_{n,i+1,1}+a_{n,i+1,2})}{(b_{n,i,1}-b_{n,i,2})},$
$\displaystyle
b_{n+1,m,2}=b_{n,0,2}+\sum_{i=1}^{n}\frac{(a_{i+1,m,1}+a_{i+1,m,2})-(a_{i,m+1,1}+a_{i,m+1,2})}{(b_{i,m,1}-b_{i,m,2})},$
(3.42)

where the right hand sides can be obtained from the moments by (3.41) and
(3.9).

###### Proof.

The relations (3.42) are obtained from the discrete zero curvature condition
(3.39) (see also (1.7)) by summation of the corresponding relations for
consecutive indices. ∎

Now we re-derive the dm-Toda equations (2.18) that we already obtained in
Subsection 3.4. However, in this case we follow the consistency approach from
[20] and [49]. In particular, we get the Lax pair here by using a method that
is the adaptation of the one from [49] (see Proposition 3.1). Since we only
consider the case of two measures, we need to consider the family of two
measures $x^{t}d\mu_{1}(x)$ and $x^{t}d\mu_{2}(x)$, where
$t\in{\mathbb{Z}}_{+}$ is the discrete time. In other words, we have two
sequences of moments $\big{\\{}s_{j}^{(1)}\big{\\}}_{j=0}^{\infty}$ and
$\big{\\{}s_{j}^{(2)}\big{\\}}_{j=0}^{\infty}$ and we consider their
truncations

$\displaystyle\big{\\{}s_{j+t}^{(1)}\big{\\}}_{j=0}^{\infty},\qquad\big{\\{}s_{j+t}^{(2)}\big{\\}}_{j=0}^{\infty},\qquad
t\in{\mathbb{Z}}_{+},$

which are actually given by the measures $x^{t}d\mu_{1}(x)$ and
$x^{t}d\mu_{2}(x)$ Clearly, these sequences of moments generate a family of
multiple orthogonal polynomials, which, as we have already seen, have the
following form

$\displaystyle
P_{n,m}^{t}(x)=\frac{1}{\tau_{n,m}^{t}}\left|\begin{matrix}\mu_{t,1}&\cdots&\mu_{t+n-1,1}\\\
\mu_{t+1,1}&\cdots&\mu_{t+n,1}\\\ \vdots&\cdots&\vdots\\\
\mu_{t+n+m,1}&\cdots&\mu_{t+2n+m-1,1}\end{matrix}\
\begin{matrix}\mu_{t,2}&\cdots&\mu_{t+m-1,2}\\\
\mu_{t+1,2}&\cdots&\mu_{t+m,2}\\\ \vdots&\cdots&\vdots\\\
\mu_{t+n+m,2}&\cdots&\mu_{t+n+2m-1,2}\end{matrix}\ \begin{matrix}1\\\ x\\\
\vdots\\\ x^{n+m}\end{matrix}\right|,$

with

$\displaystyle\tau_{n,m}^{t}=\left|\begin{matrix}\mu_{t,1}&\cdots&\mu_{t+n-1,1}\\\
\mu_{t+1,1}&\cdots&\mu_{t+n,1}\\\ \vdots&\cdots&\vdots\\\
\mu_{t+n+m-1,1}&\cdots&\mu_{t+2n+m-2,1}\end{matrix}\
\begin{matrix}\mu_{t,2}&\cdots&\mu_{t+m-1,2}\\\
\mu_{t+1,2}&\cdots&\mu_{t+m,2}\\\ \vdots&\cdots&\vdots\\\
\mu_{t+n+m-1,2}&\cdots&\mu_{t+n+2m-2,2}\end{matrix}\right|.$

As a matter of fact, we obtained an analogue of the Christoffel transformation
in the case of multiple orthogonal polynomials in Proposition 3.2.
Nevertheless, let’s do it again but this time we apply the following two
different forms of the Sylvester identity

$\displaystyle|A||A_{1,n+m+1;n+m,n+m+1}|=|A_{1;n+m}||A_{n+m+1;n+m+1}|-|A_{n+m+1;n+m}||A_{1;n+m+1}|,$
$\displaystyle|A||A_{1,n+m+1;n,n+m+1}|=|A_{1;n}||A_{n+m+1;n+m+1}|-|A_{n+m+1;n}||A_{1;n+m+1}|,$

to the determinant $\tau_{n,m}^{t}P_{n,m}^{t}(x)$. Evidently, this leads to
the relations

$\displaystyle
P_{n,m}^{t}(x)=xP_{n,m-1}^{t+1}(x)-A_{n,m-1,2}^{t}P_{n,m-1}^{t}(x),$
$\displaystyle
P_{n,m}^{t}(x)=xP_{n-1,m}^{t+1}(x)-A_{n-1,m,1}^{t}P_{n-1,m}^{t}(x),$ (3.43)

where the coefficients are defined by the formulas

$\displaystyle
A_{n,m-1,2}^{t}=\frac{\tau_{n,m-1}^{t}\tau_{n,m}^{t+1}}{\tau_{n,m}^{t}\tau_{n,m-1}^{t+1}},\qquad
A_{n-1,m,1}^{t}=\frac{\tau_{n-1,m}^{t}\tau_{n,m}^{t+1}}{\tau_{n,m}^{t}\tau_{n-1,m}^{t+1}}.$
(3.44)

Now, based on the relations (3.38) and (3.43), we can extend Proposition 3.1
to the context of multiple orthogonal polynomials. Thus, we are in the
position to complete the associated discrete integrable system (3.39) on
${\mathbb{Z}}_{+}^{2}$ to a discrete integrable system on
${\mathbb{Z}}_{+}^{3}$. To this end, we first obtain the following relations

$\displaystyle
P_{n-1,m}^{t+1}(x)=\frac{1}{x}P_{n,m}^{t}(x)+\frac{A_{n-1,m,1}^{t}}{x}P_{n-1,m}^{t}(x),$
$\displaystyle
P_{n,m-1}^{t+1}(x)=\frac{1}{x}P_{n,m}^{t}(x)+\frac{A_{n,m-1,2}^{t}}{x}P_{n,m-1}^{t}(x),$
$\displaystyle
P_{n,m}^{t+1}(x)=\left(1-\frac{b_{n,m,2}^{t}-A_{n,m,2}^{t}}{x}\right)P_{n,m}^{t}-\frac{a_{n,m,1}^{t}}{x}P_{n-1,m}^{t}-\frac{a_{n,m,2}^{t}}{x}P_{n,m-1}^{t},$
(3.45)

by manipulations with (3.43) and (3.38).

Now we see that we have a lot of options to travel over ${\mathbb{Z}}_{+}^{3}$
using the above-given relations. It is obvious that we don’t have to use all
of them to do that. However, applying different formulas when moving along the
same path leads to consistency relations and the following statement contains
all of them.

######  Theorem 3.12.

Let us consider the vector-valued wave function

$\displaystyle\Psi_{n,m,t}(x)=\big{(}P_{n,m}^{t}(x),P_{n-1,m}^{t}(x),P_{n,m-1}^{t}(x)\big{)}^{\top}.$

Then three families of matrices given by the formulas

$\displaystyle
L_{n,m,t}=\begin{pmatrix}x-b^{t}_{n,m,1}&-{a^{t}_{n,m,1}}&-{a^{t}_{n,m,2}}\\\
1&0&0\\\ 1&0&{b^{t}_{n,m-1,2}}-b^{t}_{n,m-1,1}\end{pmatrix},$ $\displaystyle
M_{n,m,t}=\begin{pmatrix}x-{b^{t}_{n,m,2}}&-{a^{t}_{n,m,1}}&-{a^{t}_{n,m,2}}\\\
1&b^{t}_{n-1,m,1}-{b^{t}_{n-1,m,2}}&0\\\ 1&0&0\end{pmatrix},$

and

$\displaystyle
N_{n,m,t}=\begin{pmatrix}1-\frac{{b^{t}_{n,m,2}}-A_{n,m,2}^{t}}{x}&-\frac{{a^{t}_{n,m,1}}}{x}&-\frac{{a^{t}_{n,m,2}}}{x}\\\
\frac{1}{x}&\frac{A_{n-1,m,1}^{t}}{x}&0\\\
\frac{1}{x}&0&\frac{A_{n,m-1,2}^{t}}{x}\end{pmatrix}$

are the transition matrices for $\Psi_{n,m,k}$ and they satisfy the following
relations

$\displaystyle 0=L_{n,m+1,t}M_{n,m,t}-M_{n+1,m,t}L_{n,m,t},$ $\displaystyle
0=M_{n,m,t+1}N_{n,m,t}-N_{n,m+1,t}M_{n,m,t},$ $\displaystyle
0=L_{n,m,t+1}N_{n,m,t}-N_{n+1,m,t}L_{n,m,t},$ (3.46)

which give the discrete zero curvature condition.

###### Proof.

To begin with, let us notice that the relations (3.38), and (3.45) can be used
to get the following vector equalities

$\displaystyle\Psi_{n+1,m,t}=L_{n,m,t}\Psi_{n,m,t},\\!\qquad\Psi_{n,m+1,t}=M_{n,m,t}\Psi_{n,m,t},\\!\qquad\Psi_{n,m,t+1}=N_{n,m,t}\Psi_{n,m,t}.\\!\\!\\!$
(3.47)

The latter system means that the matrices $L_{n,m,t}$, $M_{n,m,t}$, and
$N_{n,m,t}$ are transition matrices for the wave function $\Psi_{n,m,t}$.
Next, one can easily see that the consistency of (3.47) leads to the following
relations

$\displaystyle
L_{n,m+1,t}M_{n,m,t}\Psi_{n,m,t}=M_{n+1,m,t}L_{n,m,t}\Psi_{n,m,t},$
$\displaystyle
M_{n,m,t+1}N_{n,m,t}\Psi_{n,m,t}=N_{n,m+1,t}M_{n,m,t}\Psi_{n,m,t},$
$\displaystyle
L_{n,m,t+1}N_{n,m,t}\Psi_{n,m,t}=N_{n+1,m,t}L_{n,m,t}\Psi_{n,m,t}.$ (3.48)

Now to get (3.46) it remains to observe that the polynomials in the vector
$\Psi_{n,m,t}$ are linearly independent whenever $(n,m)$ is a normal index at
the moment $t$. Indeed, if there are numbers $\alpha_{1}$, $\alpha_{2}$, and
$\alpha_{3}$ such that

$\displaystyle\alpha_{1}P_{n,m}^{t}+\alpha_{2}P_{n-1,m}^{t}+\alpha_{3}P_{n,m-1}^{t}=0,$

then it follows by comparing the leading coefficients that $\alpha_{1}=0$.
Next, if one multiplies that relation by $x^{n-1}$ and integrates the
resulting relation with respect to $x^{k}d\mu_{1}(x)$ then one gets

$\displaystyle\alpha_{2}\int x^{n-1}P_{n-1,m}^{t}(x)x^{t}d\mu_{1}(x)=0$

since $P_{n,m-1}^{t}$ is orthogonal to $x^{n-1}$ with respect to
$x^{k}d\mu_{1}(x)$ by definition. Hence, $\alpha_{2}=0$ because of the
normality of the index $(n,m,t)$. In other words, normality means that the
determinant $\tau_{n,m}^{t}$ is nonvanishing. On the other hand, it is not so
hard to see that

$\displaystyle\int
x^{n-1}P_{n-1,m}^{t}(x)x^{t}d\mu_{1}(x)=\epsilon\tau_{n,m}^{t},$

where $\epsilon=\pm 1$. Analogously, one can get that $\alpha_{3}=0$. The
relations (3.48) reduce to (3.46), which is a discrete integrable system on
${\mathbb{Z}}_{+}^{3}$. ∎

######  Remark 3.13.

The statement of the above theorem shows that the system of difference
equations obtained from (3.46) is an integrable system in the sense of [20].
However, one might still wonder about the relation between (3.46) and (1.13).
Basically, (3.46) is a representation of (1.13) by means of a certain Lax
pair.

Indeed, this is the case and we are going to show how one can get
representative equations of (1.13) from the Lax pair representation (3.46). To
this end, let us consider entry $(1,3)$ of the relation

$\displaystyle 0=M_{n,m,t+1}N_{n,m,t}-N_{n,m+1,t}M_{n,m,t}.$

More precisely, the entry in question gives

$\displaystyle
0=-\big{(}x-b^{t+1}_{n,m,2}\big{)}\frac{{a^{t}_{n,m,2}}}{x}-{a_{n,m,2}}^{t+1}\frac{A_{n,m-1,2}^{t}}{x}+\left(1-\frac{b^{t}_{n,m+1,2}-A_{n,m+1,2}^{t}}{x}\right){a^{t}_{n,m,2}},$

which is equivalent to

$\displaystyle
b^{t+1}_{n,m,2}-b^{t}_{n,m+1,2}+A_{n,m+1,2}^{t}=\frac{a^{t+1}_{n,m,2}}{a^{t}_{n,m,2}}A_{n,m-1,2}^{t}.$

Now, taking into account the first relation in (3.44) and the second one in
(3.41), one can see that

$\displaystyle\frac{a^{t+1}_{n,m,2}}{a^{t}_{n,m,2}}A_{n,m-1,2}^{t}=A_{n,m,2}^{t},$

and, therefore, we arrive at

$\displaystyle b^{t+1}_{n,m,2}-b^{t}_{n,m+1,2}=A_{n,m,2}^{t}-A_{n,m+1,2}^{t}.$
(3.49)

Next, we know from (3.24) that

$\displaystyle b^{t}_{n,m,2}=A_{n,m,2}^{t}+\sum_{k=1}^{2}B_{n,m,k}^{t}.$

As a consequence, it follows from (3.49) that

$\displaystyle
A_{n,m,2}^{t+1}+\sum_{k=1}^{2}B_{n,m,k}^{t+1}=A_{n,m,2}^{t}+\sum_{k=1}^{2}B_{n,m+1,k}^{t},$

which is one of the equations from (1.13). To get another equation in (1.13),
let us take a look at entry $(1,2)$ of the relation

$\displaystyle 0=M_{n,m,t+1}N_{n,m,t}-N_{n,m+1,t}M_{n,m,t}.$

The entry in question gives

$\displaystyle
0=-\big{(}x-b^{t+1}_{n,m,2}\big{)}\frac{a^{t}_{n,m,1}}{x}-a^{t+1}_{n,m,1}\frac{A_{n-1,m,1}^{t}}{x}+\left(1-\frac{b^{t}_{n,m+1,2}-A_{n,m+1,2}^{t}}{x}\right)a^{t}_{n,m,1}$
$\displaystyle\hphantom{0=}{}-\frac{a^{t}_{n,m+1,1}}{x}\big{(}b^{t}_{n-1,m,1}-b^{t}_{n-1,m,2}\big{)},$

which is equivalent to

$\displaystyle
b^{t+1}_{n,m,2}-b^{t}_{n,m+1,2}+A_{n,m+1,2}^{t}=\frac{a^{t+1}_{n,m,1}}{a^{t}_{n,m,1}}A_{n-1,m,1}^{t}-\frac{a^{t}_{n,m+1,1}}{a^{t}_{n,m,1}}\big{(}b^{t}_{n-1,m,1}-b^{t}_{n-1,m,2}\big{)}.$

Now, the second relation in (1.7) gives

$\displaystyle\frac{a^{t}_{n,m+1,1}}{a^{t}_{n,m,1}}\big{(}b^{t}_{n-1,m,1}-b^{t}_{n-1,m,2}\big{)}=b^{t}_{n,m,1}-b^{t}_{n,m,2}.$

Hence, we arrive at

$\displaystyle
b^{t+1}_{n,m,2}-b^{t}_{n,m+1,2}+A_{n,m+1,2}^{t}=\frac{a^{t+1}_{n,m,1}}{a^{t}_{n,m,1}}A_{n-1,m,1}^{t}+b^{t}_{n,m,1}-b^{t}_{n,m,2}.$
(3.50)

Since by (3.24) we have

$\displaystyle b^{t}_{n,m,1}=A_{n,m,1}^{t}+\sum_{k=1}^{2}B_{n,m,k}^{t},\qquad
b^{t}_{n,m,2}=A_{n,m,2}^{t}+\sum_{k=1}^{2}B_{n,m,k}^{t},$

due to (3.49) the relation (3.50) reduces to

$\displaystyle
A_{n,m,2}^{t}-A_{n,m+1,2}^{t}+A_{n,m+1,2}^{t}=\frac{a^{t+1}_{n,m,1}}{a^{t}_{n,m,1}}A_{n-1,m,1}^{t}-A_{n,m,1}^{t}+A_{n,m,2}^{t}$

or, equivalently,

$\displaystyle{a^{t+1}_{n,m,1}}=\frac{A_{n,m,1}^{t}}{A_{n-1,m,1}^{t}}a^{t}_{n,m,1}.$
(3.51)

Next, according to (3.24) we have

$\displaystyle a^{t}_{n,m,1}=A_{n-1,m,1}^{t}B_{n,m,1}^{t}.$

Hence (3.51) reduces to

$\displaystyle A_{n-1,m,1}^{t+1}B_{n,m,1}^{t+1}=A_{n,m,1}^{t}B_{n,m,1}^{t},$

which is clearly another one from (1.13). In other words, we have reached the
following conclusion.

######  Proposition 3.14.

The second and third relations in (3.46) are one of the Lax representations of
(1.13).

###  3.5 The explicit solution of the dm-Toda equation

To illustrate our approach and provide the reader with an explicit example
when the scheme can be applied, let us recall that multiple Laguerre
polynomials of the second kind are given by the orthogonality relations

$\displaystyle\int_{0}^{\infty}x^{k}L_{n,m}^{\alpha}(x)x^{\alpha}e^{-c_{j}x}dx=0,\qquad
k=0,1,\ldots,n_{j}-1,$

for $j=1,2$, where $\alpha>-1$, $c_{1}$, $c_{2}>0$ and $c_{1}\neq c_{2}$.
Evidently, putting $t=\alpha$ and denoting $P_{n,m}^{t}(x)=L_{n,m}^{t}(x)$ we
get the polynomials with the discrete-time dynamics

$\displaystyle d\mu_{1}(t,x)=x^{t}e^{-c_{1}x}dx,\qquad
d\mu_{1}(t,x)=x^{t}e^{-c_{2}x}dx,$

and the corresponding coefficients

$\displaystyle a^{t}_{n,m,1}=\frac{n+m+t}{c_{1}^{2}}m,\qquad
a^{t}_{n,m,2}=\frac{n+m+t}{c_{2}^{2}}n,$ $\displaystyle
b^{t}_{n,m,1}=\frac{n+m+t}{c_{1}}+\frac{n}{c_{1}}+\frac{m}{c_{2}},\qquad
b^{t}_{n,m,2}=\frac{n+m+t}{c_{2}}+\frac{n}{c_{1}}+\frac{m}{c_{2}}$

are a solution of (3.46). The multiple Laguerre polynomials
$L_{n_{1},n_{2}}^{t}$ of the second kind can be obtained using the Rodrigues
formula

$\displaystyle(-1)^{n_{1}+n_{2}}\left(\prod_{j=1}^{2}c_{j}^{n_{j}}\right)x^{t}L_{n_{1},n_{2}}^{t}(x)=\prod_{j=1}^{2}\left(e^{c_{j}x}\frac{d^{n_{j}}}{dx^{n_{j}}}e^{-c_{j}x}\right)x^{n_{1}+n_{2}+t},$

where the differential operators in the product can be taken in any order
[12].

### Acknowledgements

A.I. Aptekarev was supported by grant RScF-14-21-00025. M. Derevyagin thanks
the hospitality of Department of Mathematics of KU Leuven, where his part of
the research was initiated while he was a postdoc there. M. Derevyagin and W.
Van Assche gratefully acknowledge the support of FWO Flanders project
G.0934.13, KU Leuven research grant OT/12/073 and the Belgian Interuniversity
Attraction Poles programme P07/18. H. Miki was supported by JSPS KAKENHI Grant
Number 15K17561. Also, M. Derevyagin and H. Miki are grateful to S. Tsujimoto,
L. Vinet, A. Zhedanov for valuable discussions and comments. Finally, all the
authors thank the anonymous referees for their careful reading of the
manuscript and for their remarks that improved the presentation of the paper.

## References

  * [1] Adler M., Horozov E., van Moerbeke P., The Pfaff lattice and skew-orthogonal polynomials, Int. Math. Res. Not. 1999 (1999), 569–588, solv-int/9903005.
  * [2] Adler M., van Moerbeke P., Generalized orthogonal polynomials, discrete KP and Riemann–Hilbert problems, Comm. Math. Phys. 207 (1999), 589–620, nlin.SI/0009002.
  * [3] Adler V.E., Discrete equations on planar graphs, J. Phys. A: Math. Gen. 34 (2001), 10453–10460.
  * [4] Álvarez-Fernández C., Fidalgo Prieto U., Mañas M., Multiple orthogonal polynomials of mixed type: Gauss–Borel factorization and the multi-component 2D Toda hierarchy, Adv. Math. 227 (2011), 1451–1525, arXiv:1004.3916.
  * [5] Angelesco A., Sur deux extensions des fractions continues algébriques, C. R. Acad. Sci. Paris 168 (1919), 262–265.
  * [6] Aptekarev A.I., Multiple orthogonal polynomials, J. Comput. Appl. Math. 99 (1998), 423–447.
  * [7] Aptekarev A.I., Spectral problems of high-order recurrences, in Spectral theory and differential equations, Amer. Math. Soc. Transl. Ser. 2, Vol. 233, Amer. Math. Soc., Providence, RI, 2014, 43–61.
  * [8] Aptekarev A.I., The Mhaskar–Saff variational principle and location of the shocks of certain hyperbolic equations, in Modern Trends in Constructive Function Theory, Contemporary Mathematics, Vol. 661, Amer. Math. Soc., Providence, RI, 2016, 167–186.
  * [9] Aptekarev A.I., Bleher P.M., Kuijlaars A.B.J., Large $n$ limit of Gaussian random matrices with external source. II, Comm. Math. Phys. 259 (2005), 367–389, math-ph/0408041.
  * [10] Aptekarev A.I., Branquinho A., Padé approximants and complex high order Toda lattices, J. Comput. Appl. Math. 155 (2003), 231–237.
  * [11] Aptekarev A.I., Branquinho A., Marcellán F., Toda-type differential equations for the recurrence coefficients of orthogonal polynomials and Freud transformation, J. Comput. Appl. Math. 78 (1997), 139–160.
  * [12] Aptekarev A.I., Branquinho A., Van Assche W., Multiple orthogonal polynomials for classical weights, Trans. Amer. Math. Soc. 355 (2003), 3887–3914.
  * [13] Aptekarev A.I., Derevyagin M., Van Assche W., On 2D discrete Schrödinger operators associated with multiple orthogonal polynomials, J. Phys. A: Math. Theor. 48 (2015), 065201, 16 pages, arXiv:1410.1332.
  * [14] Aptekarev A.I., Derevyagin M., Van Assche W., Discrete integrable systems generated by Hermite–Padé approximants, Nonlinearity 29 (2016), 1487–1506, arXiv:1409.4053.
  * [15] Barrios Rolanía D., Branquinho A., Foulquié Moreno A., On the relation between the full Kostant–Toda lattice and multiple orthogonal polynomials, J. Math. Anal. Appl. 377 (2011), 228–238, arXiv:0911.2856.
  * [16] Ben Cheikh Y., Douak K., On the classical $d$-orthogonal polynomials defined by certain generating functions. I, Bull. Belg. Math. Soc. Simon Stevin 7 (2000), 107–124.
  * [17] Ben Cheikh Y., Douak K., On the classical $d$-orthogonal polynomials defined by certain generating functions. II, Bull. Belg. Math. Soc. Simon Stevin 8 (2001), 591–605.
  * [18] Bernstein L., The Jacobi–Perron algorithm – Its theory and application, Lecture Notes in Math., Vol. 207, Springer-Verlag, Berlin – New York, 1971.
  * [19] Bleher P.M., Kuijlaars A.B.J., Random matrices with external source and multiple orthogonal polynomials, Int. Math. Res. Not. 2004 (2004), 109–129, math-ph/0307055.
  * [20] Bobenko A.I., Suris Yu.B., Integrable systems on quad-graphs, Int. Math. Res. Not. 2002 (2002), 573–611, nlin.SI/0110004.
  * [21] Bueno M.I., Marcellán F., Darboux transformation and perturbation of linear functionals, Linear Algebra Appl. 384 (2004), 215–242.
  * [22] Coussement J., Kuijlaars A.B.J., Van Assche W., Direct and inverse spectral transform for the relativistic Toda lattice and the connection with Laurent orthogonal polynomials, Inverse Problems 18 (2002), 923–942, math.CA/0204155.
  * [23] Coussement J., Van Assche W., An extension of the Toda lattice: a direct and inverse spectral transform connected with orthogonal rational functions, Inverse Problems 20 (2004), 297–318.
  * [24] Coussement J., Van Assche W., A continuum limit of the relativistic Toda lattice: asymptotic theory of discrete Laurent orthogonal polynomials with varying recurrence coefficients, J. Phys. A: Math. Gen. 38 (2005), 3337–3366.
  * [25] Douak K., Maroni P., Une caractérisation des polynômes $d$-orthogonaux “classiques”, J. Approx. Theory 82 (1995), 177–204.
  * [26] Fidalgo Prieto U., López Lagomasino G., Nikishin systems are perfect, Constr. Approx. 34 (2011), 297–356, arXiv:1001.0554.
  * [27] Filipuk G., Haneczok M., Van Assche W., Computing recurrence coefficients of multiple orthogonal polynomials, Numer. Algorithms 70 (2015), 519–543, arXiv:1406.0364.
  * [28] Gragg W.B., The Padé table and its relation to certain algorithms of numerical analysis, SIAM Rev. 14 (1972), 1–16.
  * [29] Grammaticos B., Kosmann-Schwarzbach Y., Tamizhmani T. (Editors), Discrete integrable systems, Lecture Notes in Phys., Vol. 644, Springer-Verlag, Berlin, 2004.
  * [30] Hermite C., Sur la fonction exponentielle, C. R. Acad. Sci. Paris 77 (1873), 18–24; 74–79; 226–233.
  * [31] Hirota R., Conserved quantities of “random-time Toda equation”, J. Phys. Soc. Japan 66 (1997), 283–284.
  * [32] Ismail M.E.H., Classical and quantum orthogonal polynomials in one variable, Encyclopedia of Mathematics and its Applications, Vol. 98, Cambridge University Press, Cambridge, 2005.
  * [33] Jack I., Jones D.R.T., Panvel J., Quantum non-abelian Toda field theories, Internat. J. Modern Phys. A 9 (1994), 3631–3656, hep-th/9308080.
  * [34] Jacobi C.G.J., Heine E., Ueber die Auslösung der Gleichung $a_{l}x_{1}+a_{2}x_{2}+\cdots+a_{n}x_{n}=f\cdot u$, J. Reine Angew. Math. 69 (1868), 1–28.
  * [35] Jacobi C.G.J., Heine E., Allgemeine Theorie der kettenbruchähnlichen Algorithmen, in welchen jede Zahl aus drei vorhergehenden gebildet wird, J. Reine Angew. Math. 69 (1868), 29–64.
  * [36] Kac M., van Moerbeke P., On an explicitly soluble system of nonlinear differential equations related to certain Toda lattices, Adv. Math. 16 (1975), 160–169.
  * [37] Mahler K., Perfect systems, Compositio Math. 19 (1968), 95–166.
  * [38] Manakov S.V., Santini P.M., Solvable vector nonlinear Riemann problems, exact implicit solutions of dispersionless PDEs and wave breaking, J. Phys. A: Math. Theor. 44 (2011), 345203, 19 pages, arXiv:1011.2619.
  * [39] Miki H., Tsujimoto S., Vinet L., Zhedanov A., An algebraic model for the multiple Meixner polynomials of the first kind, J. Phys. A: Math. Theor. 45 (2012), 325205, 11 pages, arXiv:1203.0357.
  * [40] Miki H., Vinet L., Zhedanov A., Non-Hermitian oscillator Hamiltonians and multiple Charlier polynomials, Phys. Lett. A 376 (2011), 65–69, arXiv:1106.5243.
  * [41] Moser J., Three integrable Hamiltonian systems connected with isospectral deformations, Adv. Math. 16 (1975), 197–220.
  * [42] Ndayiragije F., Van Assche W., Multiple Meixner polynomials and non-Hermitian oscillator Hamiltonians, J. Phys. A: Math. Gen. 46 (2013), 505201, 17 pages, arXiv:1310.0982.
  * [43] Nikishin E.M., Simultaneous Padé approximants, Math. USSR Sb. 41 (1982), 409–425.
  * [44] Nikishin E.M., Sorokin V.N., Rational approximations and orthogonality, Translations of Mathematical Monographs, Vol. 92, Amer. Math. Soc., Providence, RI, 1991.
  * [45] Nuttall J., Asymptotics of diagonal Hermite–Padé polynomials, J. Approx. Theory 42 (1984), 299–386.
  * [46] Papageorgiou V., Grammaticos B., Ramani A., Orthogonal polynomial approach to discrete Lax pairs for initial-boundary value problems of the QD algorithm, Lett. Math. Phys. 34 (1995), 91–101.
  * [47] Perron O., Grundlagen für eine Theorie des Jacobischen Kettenbruchalgorithmus, Math. Ann. 64 (1907), 1–76.
  * [48] Santini P.M., Nieszporski M., Doliwa A., Integrable generalization of the Toda law to the square lattice, Phys. Rev. E 70 (2004), 056615, 6 pages, nlin.SI/0409050.
  * [49] Spicer P.E., Nijhoff F.W., van der Kamp P.H., Higher analogues of the discrete-time Toda equation and the quotient-difference algorithm, Nonlinearity 24 (2011), 2229–2263, arXiv:1005.0482.
  * [50] Spiridonov V., Zhedanov A., Discrete Darboux transformations, the discrete-time Toda lattice, and the Askey–Wilson polynomials, Methods Appl. Anal. 2 (1995), 369–398.
  * [51] Spiridonov V., Zhedanov A., Spectral transformation chains and some new biorthogonal rational functions, Comm. Math. Phys. 210 (2000), 49–83.
  * [52] Suris Yu.B., The problem of integrable discretization: Hamiltonian approach, Progress in Mathematics, Vol. 219, Birkhäuser Verlag, Basel, 2003.
  * [53] Toda M., Vibration of a chain with nonlinear interaction, J. Phys. Soc. Japan 22 (1967), 431–436.
  * [54] Toda M., Theory of nonlinear lattices, Springer Series in Solid-State Sciences, Vol. 20, 2nd ed., Springer-Verlag, Berlin, 1989.
  * [55] Van Assche W., Nearest neighbor recurrence relations for multiple orthogonal polynomials, J. Approx. Theory 163 (2011), 1427–1448, arXiv:1104.3778.
  * [56] Van Iseghem J., Vector Stieltjes continued fraction and vector QD algorithm, Numer. Algorithms 33 (2003), 485–498.