Multidimensional Toda Lattices: Continuous and Discrete Time

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.


Introduction
The continuous-time Toda lattice [53,54] a n (t) = a n (t)(b n−1 (t) − b n (t)), b n (t) = a n (t) − a n−1 (t), a n > 0, t ∈ R + , n ∈ Z + , (1.1) and the discrete-time Toda lattice [31,50] A t+1 n + B t+1 n = A t n + B t n+1 , A t+1 n−1 B t+1 n = A t n B t n , t, n ∈ 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 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 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 dµ(x, t) = e −xt dµ(x), (1.3) and the functions a n (t) and b n (t) appear as the coefficients of the three-term recurrence relation xP n (x; t) = P n+1 (x; t) + b n (t)P n (x; t) + a n (t)P n−1 (x; t), n ∈ Z + , (1.4) where P n (x; t) are monic polynomials in the variable x, orthogonal with respect to dµ(x, t). Thus, the direct and inverse spectral transformations {a n (t), b n (t)} dµ(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 n , n ∈ Z r + is defined as a set of polynomials of degree | n| = n 1 + · · · + n r , satisfying the orthogonality relations L j [x i P n (x)] = x i P n (x)dµ j (x) = 0, i = 0, . . . , n j − 1, 1 ≤ j ≤ r. (1.5) This definition gives | n| homogeneous and linear equations for the | n| + 1 coefficients of the polynomial P n . In fact, we can choose P n to be monic and, so, it leads to a linear system of | 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 n is called normal. Hence, for normal indices the monic polynomial P n can be determined. If all multi-indices n ∈ Z r + on the lattice are normal then the system of measures {µ j } r j=1 (or functionals {L j } r j=1 ) generating the lattice of polynomials {P n } in (1.5) is called a perfect system. In other words, for a perfect system the polynomial P n is well defined for any n ∈ 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): xP n (x) = P n+ e i (x) + b n,i P n (x) + r k=1 a n,k P n− e k (x), i = 1, . . . , r, n ∈ Z r + , (1.6) where e i = (0, . . . , 0, 1, 0, . . . , 0) is the i-th vector of the standard basis in Z r + . Unlike in the ordinary case, the coefficients of the recurrence relations {a n,i , b n,i } for m-OPs are not independent but are required to satisfy the difference equations [55] r k=1 (a n+ e j ,k − a n+ e i ,k ) = (b n,j − b n,i )(b n+ e j ,i − b n,i ), a n,i a n+ e j ,i = b n− e i ,j − b n− e i ,i b n,j − b n,i , 1 ≤ i = j ≤ 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 dµ j (x, t) = e −tx dµ j (x, 0), 1 ≤ j ≤ r, (1.8) generates m-OPs with normal indices n and { n ± e k } r k=1 in a neighborhood of t = 0. Then the recurrence coefficients from (1.6) satisfy locally the equationṡ a n,k = a n,k [b n− e k ,k − b n,k ], b n,k = r j=1 (a n,j − a n+ e k ,j ), 1 ≤ k ≤ r. (1.9) If the system {dµ j (x, t)} r j=1 in (1.8) is perfect for t 0, then the Cauchy problem for (1.9) has a global solution which can be obtained by the direct and inverse spectral transformations {a n,k (t), b n,k (t)} {dµ 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µ j } r j=1 consists of the spectral measures of the r marginal onedimensional difference operators defined by the three-term recurrence relations (1.4) with coefficients a (j) n := a ne j ,j , b (j) n := b ne j ,j , j = 1, . . . , 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 n,k , b n,k }. Therefore, the direct and inverse spectral trans- 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 rorthogonal polynomials [16,17,25]) that is generated in the following manner 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] xq (1.11) If the system of measures {dµ j } r j=1 evolves subject to (1.8), then (see [2]) the equations for the coefficients from (1.11) can be written in the Lax pair forṁ 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, {β N } on the main diagonal and {α (k) N } r k=1 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 generates the m-OPs {P t n (x)} with normal indices n and { n ± e k } r k=1 . Then the following nonlinear dif ference system on a semi-infinite lattice can be solved locally with respect to the space variable. Moreover, if the system of measures is perfect for any discrete time t ∈ Z + , then the solution exists globally and it is given by the formulas .
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

Preliminaries
Let us briefly go over the basic definitions, see [32,Chapter 23] for details. The m-OPs P n defined in (1.5) are called multiple orthogonal (or Hermite-Padé) polynomials of type II. In what follows, we suppose that P 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 n,1 , . . . , C n,r ) are such that C n,j is a polynomial of degree at most n j −1 with the orthogonality relations r j=1 x k C n,j (x)dµ j (x) = 0, k = 0, 1, . . . , | n| − 2, and the normalization We assume that the r measures µ 1 , . . . , µ r are all absolutely continuous with respect to a measure µ and that dµ j (x) = w j (x)dµ(x) and we will use the following notation 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 n generated by type I multiple orthogonal polynomials satisfy (2.1) Furthermore, for the type II polynomials we have (1.6): xP n (x) = P n+ e k (x) + b n,k P n (x) + r j=1 a n,j P n− e j (x), (2.2) and using the type I polynomials we get .
For type I we have a similar recurrence relation xQ n (x) = Q n− e k (x) + b n− e k ,k Q n (x) + r j=1 a n,j Q n+ e j (x). (2.4) We need to point out two useful relations here. If we multiply (2.2) by Q n (x) and integrate, then the biorthogonality gives xP n (x)Q n (x)dµ(x) = r j=1 a n,j . (2.5) The orthogonality properties of P n imply that P n (x)Q n+ e k (x)dµ(x) = P n C n+ e k ,k (x)dµ k (x), so we arrive at 1 = κ n+ e k ,k P n (x)x n k dµ k (x). (2.6)

Time dynamics
In this subsection we present a proof of Theorem 1.1. To this end, we consider a perfect system of measures (µ 1 (x, t), . . . , µ r (x, t)) that depend on time as in (1.8), i.e., dµ k (x, t) := e −tx dµ k (x) for 1 ≤ k ≤ r. Now the corresponding type I and type II multiple orthogonal polynomials depend on time as well, that is, we have Q n (x; t) and P n (x; t). Our goal is to show that the relations (2.3) evaluated at the moment t b n,k (t) = xP n (x; t)Q n+ e k (x; t)e −xt dµ(x), a n, imply multiple Toda (m-Toda) equations (1.9), i.e., for 1 ≤ k ≤ ṙ a n,k = a n, andḃ n,k = r j=1 (a n,j − a n+ e k ,j ). (2.8) From the biorthogonality (2.1) we get P n (x; t)Q n+ e k (x; t)e −xt dµ(x) = 1.
If we take derivatives with respect to t, then this gives Ṗ n (x; t)Q n+ e k (x; t)e −xt dµ(x) + P n (x; t)Q n+ e k (x; t)e −xt dµ(x) Observe thatṖ n is a polynomial of degree at most | n| − 1 since P n is a monic polynomial, hence the orthogonality for type I multiple orthogonal polynomials gives Ṗ n (x; t)Q n+ e k (x; t)e −xt dµ(x) = 0.
The orthogonality for type II multiple orthogonal polynomials gives P n (x; t)Q n+ e k (x; t)e −xt dµ(x) = P n (x; t)Ȧ n+ e k ,k (x; t)e −xt dµ k (x). (2.9) If we use C n+ e k ,k (x; t) = κ n+ e k ,k (t)x n k + · · · , then we find P n (x; t)Ȧ n+ e k ,k (x; t)e −xt dµ k (x) = κ n+ e k ,k x n k P n (x; t)e −xt dµ k (x) =κ where the last equality follows from (2.5 κ n+ e k ,k =ȧ n,k a n,k , a n+ 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 b n,k (t) = xP n (x; t)Q n+ e k (x; t)e −xt dµ(x).
Taking the derivative with respect to t giveṡ b n,k (t) = xṖ n (x; t)Q n+ e k (x; t)e −xt dµ(x) If we use (2.4) then we find xṖ n (x; t)Q n+ e k (x; t)e −xt dµ(x) a n+ e k ,j Q n+ e j + e k (x; t)   e −xt dµ(x).
SinceṖ n is of degree at most | n|−1, the orthogonality for type I multiple orthogonal polynomials gives Taking the derivative of hence xṖ n (x; t)Q n+ e k (x; t)e −xt dµ(x) = r j=1 a n,j (t).
If we use (2.2) then we find The orthogonality of type II multiple orthogonal polynomials gives From (2.9) and (2.10) we recall that If we take the derivative of and after using (2.2) and the biorthogonality, we find This gives Finally use (2.2), (2.4) and the biorthogonality to find a n+ e k ,j + b 2 n,k + r j=1 a n,j .
Combining all these results then giveṡ b n,k = r j=1 a n,j − r j=1 a n+ e k ,j , which is the same as (2.8). This finishes the proof of Theorem 1.1.

Time dependent m-OPs
Introducing the moments defined by the functional (1.5) we have a determinant expression for m-OPs From this expression, it is easy to find that the multi-index n ∈ Z n + is normal iff τ n = 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 n,j , b n,j } is directly verified 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 This transformation can also be interpreted in terms of the linear functionals (1.5): 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 {µ i,j } with (2.14)) are given, then the following relation for the corresponding m-OPs holds d dt P n (x) = r k=1 a n,k P n− e k (x), (2.16) where {a n,j } are the coefficients of the recurrence relation in (1.6).
Proof . To begin with, introduce the notation τ n,x := 0 n 1 , 1 n 1 , . . . , (n 1 − 1) n 1 , . . . , 0 nr , . . . , (n r − 1) nr , x Then τ n and P n (x) can be rewritten as where e = (0, . . . , 0, 1) T . From the relation (2.14), it is easy to find that τ n and τ n,x are Wronskian matrices, which amounts to Using these notations and relations, we can calculate the derivative of P n (x): In the calculation of (2.17), we have used the Plücker relation, a well-known identity for determinants 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") xP n (x) = P n+ e j (x) + b n,j P n (x) + r k=1 a n,k P n− e k (x), j = 1, . . . , r, d dt P n (x) = r k=1 a n,k P n− 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: Recalling the determinant expression of the solution (2.13), one can easily find that the coefficients a n,j and b n,j are also expressed as follows Substituting this expression into (2.18), we can can get the following statement.
Proposition 2.2. The τ -function verif ies the following bilinear equation This bilinear equation is also a generalization of that of the ordinary Toda equation.

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 L 1 , 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 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), 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 with the contiguous relations for the initial values We get the 2-orthogonal polynomials {q n } in the following manner It can easily be checked that {q n (x)} satisfy the four-term recurrence relation where the coefficients are 2n = a n,n,1 + a n,n,2 , α Taking all this into account, it is straightforward to find the spectral transformation of the 2-orthogonal polynomials q n as followṡ 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 is satisfied. Moreover, one can rewrite the system in the Lax form 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 ≥ 2.

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 where δ > −1, κ 1 = κ 2 and t > −κ 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 b m,n,1 = 2m + n + δ + 1 It is easily verified that the corresponding linear functionals L 1 , 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 Interestingly, it is obvious that these solutions have a pole at t = −κ 1 and t = −κ 2 , which shows the presence of singularities at certain finite times.
3 The discrete-time higher analogues of the Toda lattice

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µ on (0, +∞) for which all the moments exist. Clearly, the measure x t dµ(x), which is defined on (0, +∞), is also positive for all t ∈ Z + . Let P t n be the family of polynomials orthogonal with respect to the measure x t dµ(x) on (0, +∞). Introducing the moments µ j of the measure dµ one can easily check that the monic orthogonal polynomials P t n can be presented in the following manner with the corresponding Hankel determinant 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 |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 × 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 α;β denotes the submatrix of A that is obtained from A by removing the αth column and βth row. Next, applying two different forms of the Sylvester identity |A||A 1,n+1;1,n+1 | = |A 1;1 ||A n+1;n+1 | − |A 1;n+1 ||A n+1;1 |, |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 τ t n P t n leads to the relations The transformation (3.1) from P t n to P t+1 n 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µ(x) provided that the polynomials orthogonal with respect to dµ(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 t+1 n can be represented by means of the Christoffel-Darboux kernel in the following manner 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]) On the other hand, if we plug (3.3) into (3.4) we get xP t+1 n (x) = P t+1 n+1 (x) + A t n + B t n+1 P t+1 n (x) + A t n B t n P t+1 n−1 .

Now, comparing the corresponding coefficients leads to (1.2):
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 Ψ n,t (x) = P t n (x), P t+1 n (x) .
Next, using (3.1) and (3.2) we derive Ψ n+1,t = L n,t Ψ n,t , Ψ n,t+1 = M n,t Ψ n,t , where the transition matrices are defined as follows The consistency of the linear systems (3.5) is then equivalent to the zero curvature condition 0 = L n,t+1 M n,t − M n+1,t L n,t , (3.6) which can be simplified to the quotient-difference scheme 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 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]) Hence, we also have formulas for the coefficients a t n and b t n in terms of Hankel determinants 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 Ψ n,t (x) = P t n (x), P t n−1 (x) .
Then the corresponding transition matrices are

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 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 it is clear that is equivalent to (3.6) and, in turn, reduces to (3.7).

Christof fel 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.
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 we can obtain the chain of m-OPs Remark 3.3. We can derive the nonlinear equations from the compatibility condition of the relations (3.15) themselves 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.
where a t n,k are the coefficients of the nearest neighbor recurrence relation Proof . First, consider the polynomial P t n (x) − P t+1 n (x), which is a polynomial of degree | n| − 1. From the multiple orthogonality relation and (3.14), we can easily check This readily shows P t n (x) − P t+1 n (x) is orthogonal to all polynomials of degree less than n j − 1 with respect to the linear functional L t+1 j . Hence, we can write P t n (x) − P t+1 n (x) as a linear combination of the polynomials P t+1 n− e j , j = 1, . . . , r, which form a basis for the linear space of all polynomials of degree less than | n| and satisfy the multiple orthogonality conditions n (x) = 0, j = 1, . . . , r, for k ≤ n j − 2. We write P t n (x) − P t+1 n (x) = r k=1 B t n,k P t+1 n− e k (x). From (3.19), we can easily find , k = 1, . . . , r. (3.20) Here, some calculations show We can also calculate the coefficients of (3.18) from the multiple orthogonality relation .
Remark 3.5. Although we can formally derive all the relations in Section 3 for λ t that changes with the discrete time t, we are only concerned with the case λ t = λ and, particularly, λ t = 0. The reason is that the basis for our construction of integrable systems is a perfect system of measures dµ 1 , . . . , dµ 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 λ t = 0, which can be easily modified to λ t = λ (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 A t n,j , B t n,j b t n,j = A t n,j + r k=1 B t n,k + λ t , a t n,j = A t n− e j ,j B t n,j , j = 1, . . . , 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 t n,j , B t n,j }. Corollary 3.6. If we put λ t = 0, the coefficients A t n,j , B t n,j in (3.15) and (3.17) satisfy the following difference equations on n A t n,i A t n+ e i ,j = A t n,j A t n+ e j ,i , for all t and i, j = 1, . . . , 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 (x − λ t )P t+1 n (x) = P t n+ e j (x) + A t n,j P t n (x), j = 1, . . . , r, P t n (x) = P t+1 n (x) + r k=1 B t n,k P t+1 n− e k (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: 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.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 τ -function τ t n defined by where µ t i,j := L t j [x i ]. From (3.14), we get the following relation From the determinant expression of the m-OPs (2.12) we can get by means of elementary transformations of determinants We shall consider the bilinear equations for dm-Toda lattice (3.27) for the case λ t = λ. From (3.16), the dependent variable A t n,j obeys a discrete KP equation, which reduces to the following Hirota-Miwa equation

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 λ t = λ. We denote the corresponding m-OPs by {P t m,n } and we write the discrete Lax set of m-OPs as follows Let us introduce the sequence of 2-orthogonal polynomials {q t n } by the correspondence q t 2n (x) = P t n,n (x), q t 2n+1 (x) = P t n+1,n (x).
In a fashion similar to the previous section, we can also directly obtain the discrete spectral transformation of 2-orthogonal polynomials as follows From the compatibility condition of (3.34) it follows that the discrete integrable system associated with 2-orthogonal polynomials is given by

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.
x n+m provided that τ 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 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), 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 ≥ 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]. 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] a n,m,1 = τ n+1,m τ n−1,m τ n,m 2 , a n,m,2 = τ n,m+1 τ n,m−1 τ n,m (a n+1,i,1 + a n+1,i,2 ) − (a n,i+1,1 + a n,i+1,2 ) 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µ 1 (x) and x t dµ 2 (x), where t ∈ Z + is the discrete time. In other words, we have two sequences of moments s which are actually given by the measures x t dµ 1 (x) and x t dµ 2 (x) Clearly, these sequences of moments generate a family of multiple orthogonal polynomials, which, as we have already seen, have the following form P t n,m (x) = 1 τ t n,m µ t, 1 · · · µ t+n−1,1 µ t+1,1 · · · µ t+n,1 . . . · · · . . . µ t+n+m,1 · · · µ t+2n+m−1,1 µ t,2 · · · µ t+m−1,2 µ t+1,2 · · · µ t+m,2 . . . · · · . . .
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 |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 |, |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 τ t n,m P t n,m (x). Evidently, this leads to the relations P t n,m (x) = xP t+1 n,m−1 (x) − A t n,m−1,2 P t n,m−1 (x), P t n,m (x) = xP t+1 n−1,m (x) − A t n−1,m,1 P t n−1,m (x), (3.43) where the coefficients are defined by the formulas 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 Z 2 + to a discrete integrable system on Z 3 + . To this end, we first obtain the following relations by manipulations with (3.43) and (3.38). Now we see that we have a lot of options to travel over 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.  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 Ψ n+1,m,t = L n,m,t Ψ n,m,t , Ψ n,m+1,t = M n,m,t Ψ n,m,t , Ψ n,m,t+1 = N n,m,t Ψ 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 Ψ n,m,t . Next, one can easily see that the consistency of (3.47) leads to the following relations L n,m+1,t M n,m,t Ψ n,m,t = M n+1,m,t L n,m,t Ψ n,m,t , M n,m,t+1 N n,m,t Ψ n,m,t = N n,m+1,t M n,m,t Ψ n,m,t , L n,m,t+1 N n,m,t Ψ n,m,t = N n+1,m,t L n,m,t Ψ n,m,t . (3.48) Now to get (3.46) it remains to observe that the polynomials in the vector Ψ n,m,t are linearly independent whenever (n, m) is a normal index at the moment t. Indeed, if there are numbers α 1 , α 2 , and α 3 such that α 1 P t n,m + α 2 P t n−1,m + α 3 P t n,m−1 = 0, then it follows by comparing the leading coefficients that α 1 = 0. Next, if one multiplies that relation by x n−1 and integrates the resulting relation with respect to x k dµ 1 (x) then one gets α 2 x n−1 P t n−1,m (x)x t dµ 1 (x) = 0 since P t n,m−1 is orthogonal to x n−1 with respect to x k dµ 1 (x) by definition. Hence, α 2 = 0 because of the normality of the index (n, m, t). In other words, normality means that the determinant τ t n,m is nonvanishing. On the other hand, it is not so hard to see that x n−1 P t n−1,m (x)x t dµ 1 (x) = τ t n,m , where = ±1. Analogously, one can get that α 3 = 0. The relations (3.48) reduce to (3.46), which is a discrete integrable system on 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 ( and, therefore, we arrive at b t+1 n,m,2 − b t n,m+1,2 = A t n,m,2 − A t n,m+1,2 . Next, according to (3.24) we have a t n,m,1 = A t n−1,m,1 B t n,m,1 .

Hence (3.51) reduces to
A t+1 n−1,m,1 B t+1 n,m,1 = A t n,m,1 B t n,m,1 , 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).

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 ∞ 0 x k L α n,m (x)x α e −c j x dx = 0, k = 0, 1, . . . , n j − 1, for j = 1, 2, where α > −1, c 1 , c 2 > 0 and c 1 = c 2 . Evidently, putting t = α and denoting P t n,m (x) = L t n,m (x) we get the polynomials with the discrete-time dynamics e c j x d n j dx n j e −c j x x n 1 +n 2 +t , where the differential operators in the product can be taken in any order [12].