Discrete integrable systems generated by Hermite-Padé approximants
aa r X i v : . [ m a t h . C A ] S e p DISCRETE INTEGRABLE SYSTEMS GENERATED BYHERMITE-PAD´E APPROXIMANTS
ALEXANDER I. APTEKAREV, MAXIM DEREVYAGIN, AND WALTER VAN ASSCHE
Abstract.
We consider Hermite-Pad´e approximants in the framework of dis-crete integrable systems defined on the lattice Z . We show that the conceptof multiple orthogonality is intimately related to the Lax representations forthe entries of the nearest neighbor recurrence relations and it thus gives rise toa discrete integrable system. We show that the converse statement is also true.More precisely, given the discrete integrable system in question there exists aperfect system of two functions, i.e., a system for which the entire table ofHermite-Pad´e approximants exists. In addition, we give a few algorithms tofind solutions of the discrete system. Introduction
Nowadays modern technologies allow us to handle an enormous amount of in-formation. As a consequence of this development, it is in many instances moreadvantageous to face the analysis of discrete data rather than continuous data. Forthis reason we are witnessing that the world in this century requires more and morethe understanding of discrete models and that is why we decided to concentrate ourattention on studying discrete models that take their origin in orthogonality , oneof the basic mathematical concepts. Mathematically speaking, the discrete modelswe are going to consider are systems of difference equations . Recent advances ina number of mathematical fields reveal that discrete systems are in many respectseven more fundamental than continuous ones (for instance see [25], [28]).In this paper we follow the streamline of discrete integrable systems (see [8], [9])and our main interest is in discrete systems on Z represented by a field of squareinvertible matrices(1.1) L n,m , M n,m ∈ C d × d , n, m ∈ Z , which satisfy the discrete zero curvature condition (or form a Lax pair ) on Z :(1.2) L n,m +1 M n,m − M n +1 ,m L n,m = 0 . The elements of the discrete system (1.1) are transition matrices which define theevolution of the wave function Ψ n,m (1.3) Ψ n +1 ,m ( z ) = L n,m ( z )Ψ n,m ( z ) , Ψ n,m +1 ( z ) = M n,m ( z )Ψ n,m ( z ) , Date : September 6, 2018.1991
Mathematics Subject Classification.
Primary 42C05, 37K10; Secondary 37K60, 39A12,39A14, 65Q10.
Key words and phrases.
Multiple orthogonal polynomials, discrete integrable system, discretezero curvature condition, ordinary and partial difference equations, two-dimensional Schur-Euclidalgorithm, two-dimensional continued fractions, recurrence relations. and the condition (1.2) describes the consistency or integrability of the equations(1.3). In turn, the relations (1.2) represent a nonlinear system of difference equa-tions.Our findings are mainly inspired by the connection between discrete integrablesystems, orthogonal polynomials , and
Pad´e approximants (see [20], [29]). For ex-ample, the discrete dynamics(1.4) x k dµ ( x ) , k ∈ Z + of the measure dµ supported on [0 , + ∞ ) generates a family of orthogonal polyno-mials { P ( k ) n ( x ) } , deg P ( k ) n = n . These polynomials also appear as numerators ofPad´e approximants in the Pad´e table and some nearest neighbour polynomials P ( k ) n are related P ( k ) n +1 ( x ) = xP ( k +1) n ( x ) − V ( k ) n P ( k ) n , P ( k ) n +1 ( x ) = xP ( k +2) n ( x ) − W ( k ) n P ( k +1) n , where the coefficients of the relations are expressed by means of the Hankel deter-minants V ( k ) n = S ( k +1) n +1 S ( k ) n S ( k +1) n S ( k ) n +1 , W ( k ) n = S ( k +1) n +1 S ( k +1) n S ( k ) n +1 S ( k +2) n , S ( k ) n +1 = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s k . . . s n + k ... . . . ... s n + k . . . s n + k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , and s j = R ∞ x j dµ ( x ). Finally, the consistency of these relations gives the discretezero curvature condition (1.2) for the discrete system (1.1) of 2 × L n,k = − V ( k ) n x − V ( k ) n x + W ( k ) n − V ( k +1) n ! , M n,k = 1 x (cid:18) x − V ( k ) n x + W ( k ) n (cid:19) . Recall that in the theory of Pad´e approximants this discrete system becomes thequotient-difference algorithm, and in integrable systems theory it leads to thediscrete-time Toda equation (see, e.g., [31]).In the present paper we introduce a new class of discrete integrable systems of3 × Hermite-Pad´e rational approximants , which were introduced by Hermite [15] inconnection to his outstanding proof of the transcendence of e . These days thistheory is known to play an important role in various fields ranging from numbertheory [4], [6], [33] to random matrix theory [19], [7].To proceed, let us briefly consider the concept of Hermite-Pad´e rational approxi-mants (for details, see the surveys [5], [32]). Let ~f = ( f , f ) be a vector of Laurentseries at infinity(1.5) f j ( z ) = ∞ X k =0 s j,k z k +1 , j = 1 , . The
Hermite-Pad´e rational approximants (of type II) π ~n = Q (1) ~n P ~n , Q (2) ~n P ~n ! for the vector ~f and multi-index ~n = ( n , n ) ∈ N are defined bydeg P ~n ≤ | ~n | = n + n , ISCRETE INTEGRABLE SYSTEMS GENERATED BY HP APPROXIMANTS 3 (1.6) f j ( z ) P ~n ( z ) − Q ( j ) ~n ( z ) =: R ( j ) ~n ( z ) = O (cid:18) z n j +1 (cid:19) , z → ∞ , where the Q ( j ) ~n are polynomials, for j = 1 ,
2. This definition is equivalent to ahomogeneous linear system of equations for the coefficients of the polynomial P n ,n .This system always has a solution, but the solution is not necessarily unique. Inthe case of uniqueness (up to a multiplicative constant) and in case any non-trivialsolution has full degree deg P n ,n = n + n , the multi-index ( n , n ) is called normal and the polynomial P ~n can be normalized to be monic.Clearly these polynomials can be put in a table { P n,m } . If all indices of thistable are normal, then the system of functions (1.5) is called a perfect system . Thenotion of perfect systems was introduced by Mahler [21]. For perfect systems, theHermite-Pad´e polynomials (1.6) satisfy a system of recurrence relations(1.7) ( P n +1 ,m ( x ) = ( x − c n,m ) P n,m ( x ) − a n,m P n − ,m ( x ) − b n,m P n,m − ( x ) ,P n,m +1 ( x ) = ( x − d n,m ) P n,m ( x ) − a n,m P n − ,m ( x ) − b n,m P n,m − ( x ) , with a ,m = b n, = 0 for all n, m ≥ a n,m = 0, b n,m = 0 for all other indices( n, m ). As we will see, the consistency of these relations also leads to Lax pairrepresentations, where the corresponding matrices L n,m and M n,m have the forms(1.8) L n,m = x + α (1) n,m α (2) n,m α (3) n,m α (4) n +1 ,m α (5) n +1 ,m , M n,m = x + β (1) n,m α (2) n,m α (3) n,m α (4) n,m +1 α (5) n,m +1 , where the entries are related to the coefficients of the recurrence relations (1.7) forthe Hermite-Pad´e polynomials of the functions f and f as follows: c n,m = − α (1) n,m , d n,m = − β (1) n,m , a n,m = − α (4) n,m α (2) n,m , b n,m = − α (5) n,m α (3) n,m . Both sets of coefficients of the relations (1.7) and of entries of the matrices (1.8) canbe represented by the power series coefficients of the perfect system of functions(1.5)(1.9) ( a, b, c, d ) n,m ←− { s j,k } j =1 , −→ (cid:0) α (1) , β (1) , α (2) , . . . , α (5) (cid:1) n,m (details will be given below). The main message of this note is to show that there arediscrete integrable systems related to Hermite-Pad´e approximation and the theoryof such approximants allows to construct solutions of these systems provided theinitial boundary data are given. In this paper we restrict ourselves to the simplestcase of Hermite-Pad´e approximants generated by two functions, for which we havethe following result. Theorem 1.1.
The zero curvature condition (1.2) holds for a family of × transition matrices L n,m and M n,m of the form (1.8) if and only if there is a perfectsystem of two functions (1.5) such that P n,m are the Hermite-Pad´e polynomials withthe coefficients of the recurrence relations (1.7) corresponding to (1.9) . The proof of this statement is immediate from Proposition 4.2, which is a slightgeneralization of the result from [34], and Proposition 4.5, which is the conversestatement and is the basis for further development of the present paper.Since the transition matrices are explicitly known, it is not so difficult to re-rewrite the discrete zero curvature condition in terms of the coefficients of (1.7).
ALEXANDER I. APTEKAREV, MAXIM DEREVYAGIN, AND WALTER VAN ASSCHE
Thus, we have arrived at the following statement, which is simply an equivalentform of Theorem 1.1.
Theorem 1.1 ′ . The discrete Lax pair equations (1.2) for the matrices L n,m , M n,m of the form (1.8) are equivalent to the nonlinear system of difference equations forthe coefficients of the recurrence relations (1.7)(1.10) c n,m +1 = c n,m + ( a + b ) n +1 ,m − ( a + b ) n,m +1 ( c − d ) n,m d n +1 ,m = d n,m + ( a + b ) n +1 ,m − ( a + b ) n,m +1 ( c − d ) n,m a n,m +1 = a n,m ( c − d ) n,m ( c − d ) n − ,m b n +1 ,m = b n,m ( c − d ) n,m ( c − d ) n,m − n, m ≥ , where the initial data ( c, a ) n, , ( d, b ) ,m are supposed to be given and the coefficientsalso satisfy the boundary conditions a ,m = 0 = b n, . Moreover, the system has asolution such that a n,m = 0 , b n,m = 0 for n, m > if and only if the initial datacorrespond to a perfect system. This theorem explicitly gives us the underlying boundary value problem. More-over, the way it is written, one can easily get an idea about the continuum limit ofthis discrete system.Once the boundary value problem is given, it is natural to look for its solution.To this end, we present here a new method of solving the discrete system usingbranched continued fractions related to Hermite-Pad´e approximants. Furthermore,these branched continued fractions are introduced here for the first time and theyappear as the multiple orthogonal adaptation of the classical continued fractionrepresentation (3.3) that solves some inverse problems for finite Jacobi matrices[14]. Also, it is used in dynamical systems and asymptotic analysis.Another natural question is whether the initial data lead to a well-posed problem.In order to understand this issue a bit more deeper, let us notice that Hermite-Pad´eapproximants are intimately related to the notion of multiple orthogonal polyno-mials. If the coefficients of the Laurent series (1.5) are the moments of positivemeasures µ and µ supported on R (1.11) f j ( z ) = ∞ X k =0 s j,k z k +1 = Z R dµ j ( x ) z − x , s j,k = Z R x k dµ j ( x ) , then the Hermite-Pad´e denominators P n ,n from (1.6) satisfy(1.12) Z P n ,n ( x ) x k dµ j ( x ) = 0 , k = 0 , , . . . , n j − , j = 1 , . Polynomials defined by the system of orthogonality relations (1.12) are called mul-tiple orthogonal polynomials.
The idea of this concept is the following: given twomeasures ( µ , µ ), we distribute the n + n orthogonality relations between thesemeasures and aim to find a monic polynomial P n ,n of degree deg P n ,n = n + n that is orthogonal to the n first monomials with respect to the first measure andto the n first monomials with respect to the second measure.The multiple orthogonal polynomials (i.e., the Hermite-Pad´e polynomials) in-herit all the remarkable properties for Hermite-Pad´e approximants, like existence ISCRETE INTEGRABLE SYSTEMS GENERATED BY HP APPROXIMANTS 5 of the monic polynomials of full degree for the normal indices and the recurrencerelations (1.7), which were obtained for the first time in [34] for the multiple or-thogonal polynomials. In the context of our paper we use these polynomials togenerate a general class of perfect systems for which the corresponding tables ofmultiple orthogonal polynomials exist entirely.One should not think that any two measures form a perfect system. It is actuallynot a trivial task to check whether one can obtain such a table for any two measures.Moreover, there are measures for which it is impossible to define polynomials for allindices. So, in order to address this issue we give a short review of Angelesco andNikishin systems at the end of Section 5, which are actually pairs of measures forwhich the corresponding table of multiple polynomials exist entirely and, therefore,the boundary data obtained from such systems lead to well-posed problems forthe discrete system in question. Some algorithms for solving the boundary valueproblem will be mentioned in Section 5 as well.To conclude the introduction we want to remark that it is also possible to consideran analogue of the dynamics (1.4) for multiple orthogonal polynomials and get adiscrete integrable system similar to the qd-algorithm. Moreover, the qd-algorithmwas already partially adapted to the case of multiple orthogonal polynomial in [17].This will be considered in more detail elsewhere.
Structure of the paper.
The following two sections serve as an introductionto the topic. In particular we give in Section 2 more explanations about generaldiscrete integrable systems. Then, in Section 3, we consider some known 2 × × × ′ . Finally, in Section5 we provide the reader with a generic class of perfect systems such as Angelescoand Nikishin systems. These systems generate the boundary data for which thediscrete integrable system is solvable. Acknowledgements.
A.I. Aptekarev was supported by grant RScF-14-21-00025.M. Derevyagin thanks the hospitality of the Department of Mathematics of KULeuven, where his part of the research was mainly done while he was a postdocthere. M. Derevyagin and W. Van Assche gratefully acknowledge the support ofFWO Flanders project G.0934.13, KU Leuven research grant OT/12/073 and theBelgian Interuniversity Attraction Poles programme P07/18.2.
The generic Lax representations
Here we recall some basic notions in the theory of discrete integrable systemsfollowing [10] (see also [1], [9], [20], and [29]).Let us consider a regular square lattice Z , that is the set of all pairs ( n, m )of integer numbers n and m . The main object of our study is wave functions Ψ n,m defined on all the vertices ( n, m ) of Z and having their values in C k × k (forsimplicity, we restrict ourselves here to the cases k = 2 and k = 3). The wavefunction Ψ n,m depends on a complex parameter z , which is interpreted as the ALEXANDER I. APTEKAREV, MAXIM DEREVYAGIN, AND WALTER VAN ASSCHE spectral parameter. We assume that for any oriented edge the values of the wavefunction at the vertices that this edge connects are related via transition matrices L n,m and M n,m as followsΨ n +1 ,m ( z ) = L n,m ( z )Ψ n,m ( z ) , Ψ n,m +1 ( z ) = M n,m ( z )Ψ n,m ( z ) . We always require that the transition matrices are invertible and therefore one hasΨ n,m ( z ) = L − n,m ( z )Ψ n +1 ,m ( z ) , Ψ n,m ( z ) = M − n,m ( z )Ψ n,m +1 ( z ) . It is clear that the value of the wave function must not depend on the path onetakes to get to the corresponding vertex. Thus, in order that the wave functionΨ n,m is well defined, the following zero curvature condition must be satisfied(2.1) L n,m +1 M n,m − M n +1 ,m L n,m = 0 , n, m ∈ Z . As is known, the zero curvature condition is equivalent to integrability. Thus, adiscrete system that admits the representation (2.1) is called integrable [1], [8], [9].Before going further, let us take a careful look at the zero curvature condition.Observe at first that one can rewrite (2.1) as L − n,m M − n +1 ,m L n,m +1 M n,m = I. Next, from the following picture(n,m) M n,m (n,m+1) L n,m +1 (n+1,m+1) M − n +1 ,m (n+1,m) L − n,m we see that the zero curvature condition implies that the product of the transitionmatrices along the oriented simple square path on Z beginning at the vertex ( n, m )is the identity matrix. This observation can be immediately extended to the caseof domino paths by reducing them to the just considered simplest case.Now it is clear how to generalize this statement to the case of any closed orientedpath on Z . Thus the condition (2.1) means that if one fixes a closed orientedpath on the lattice Z , then the product of the transition matrices in the orderthey appear along the path must be the identity matrix. Note that this property ISCRETE INTEGRABLE SYSTEMS GENERATED BY HP APPROXIMANTS 7 resembles the Cauchy theorem for holomorphic functions and, therefore, it can beconsidered as its noncommutative multiplicative analogue for functions on Z . Therelation (2.1) is also called the Lax representation and this is one of the ways tosay that the underlying discrete system is integrable.It turns out that different types of wave functions appear in the theory of or-thogonal polynomials and they are very useful to achieve a big variety of goals.However, it has to be pointed out that the discrete integrable systems and wavefunctions, that have their origin in orthogonality, are naturally defined on N , where N = { , , , . . . } . Luckily one can appropriately extend them to Z or even to Z depending on the needs.3. Orthogonal polynomials via × matrix polynomials In this section we briefly go over the ideas in the theory of ordinary orthogonalpolynomials in order to consider more general objects in Section 4 and to get adiscrete integrable system whose Lax pair is expressed via 3 × The Schur-Euclid algorithm.
Suppose we are given a nontrivial Borel mea-sure dµ on the real line R . Assume also that all the moments of the measure dµ are finite. Then the Schur algorithm, which is a straightforward generalization ofEuclid’s algorithm, leads to the following continued fraction ϕ ( z ) = Z R dµ ( t ) t − z ∼ − z − a − b z − a − b . . . , where b j > a j ∈ R for j = 0 , , . . . . This continued fraction is called a J -fraction and such a representation actually exist for a larger class of analyticfunctions.It is natural to consider continued fractions as infinite sequences of linear frac-tional transformations. In particular, in the case of the J -fraction, one has thefollowing sequence ϕ j ( z ) = T j ( ϕ j +1 ( z )) = − z − a j + b j ϕ j +1 ( z ) , j ∈ Z + , with the initial condition ϕ = ϕ . Also, it is well known that a linear fractionaltransformation can be represented as a 2 × T j
7→ W j ( z ) = − b j b j z − a j b j ! , j ∈ Z + . Let us now introduce matrices corresponding to the approximants for the J -fraction,that is, the finite truncations of the continued fraction:(3.1) W [ n, ( z ) = W ( z ) W ( z ) . . . W n ( z ) , n ∈ Z + . Before showing how to construct a set of a wave functions and transition matriceson Z let us see what the elements of the matrix polynomial W [ n, are. To thisend, we put (cid:18) − Q P (cid:19) := (cid:18) (cid:19) , (cid:18) − Q j +1 ( z ) P j +1 ( z ) (cid:19) := W [ j, ( z ) (cid:18) (cid:19) , j ∈ Z + . ALEXANDER I. APTEKAREV, MAXIM DEREVYAGIN, AND WALTER VAN ASSCHE
Then taking into account the relation W [ j, ( z ) = W [ j − , ( z ) W j ( z ) we also havethat W [ j, ( z ) (cid:18) (cid:19) = W [ j − , ( z ) (cid:18) b j (cid:19) = (cid:18) − b j Q j ( z ) b j P j ( z ) (cid:19) , j ∈ N . So, the matrix W [ j, has the following form W [ j, ( z ) = (cid:18) − b j Q j ( z ) − Q j +1 ( z ) b j P j ( z ) P j +1 ( z ) (cid:19) , j ∈ Z + . Furthermore, rewriting the relation entrywise (cid:18) − Q j +1 ( z ) P j +1 ( z ) (cid:19) = W [ j − , ( z ) W j ( z ) (cid:18) (cid:19) = 1 b j W [ j − , ( z ) (cid:18) − z − a j (cid:19) , j ∈ N , we see that the polynomials P j , Q j are solutions of the following three-term recur-rence relation(3.2) b j − u j − ( z ) + a j u j ( z ) + b j u j +1 ( z ) = zu j ( z ) , j ∈ N , with the initial conditions P ( z ) = 1 , P ( z ) = z − a b ,Q ( z ) = 0 , Q ( z ) = 1 b . Thus the entries of the matrix W [ n, are orthogonal polynomials of the first andsecond kind and the corresponding orthogonality measure is µ . It is worth men-tioning that such 2 × − b j P j +1 ( z ) P j ( z ) = − z − a j − b j − z − a j − − b j − . . . − b z − a that will be generalized to the case of multiple orthogonal polynomials and thenit will be used in the scheme of finding solutions of the discrete system underconsideration.3.2. Riemann-Hilbert problems.
Here we consider a different interpretation ofthe Schur-Euclid algorithm in the context of Riemann-Hilbert problems, whichturned out to be quite efficient for asymptotic analysis. Recall that in [13] a fas-cinating characterization of orthogonal polynomial in terms of a Riemann-Hilbertproblem was found. We will explain this characterization here briefly. To this end,we consider a weight function w on R that is smooth and has sufficient decay at ±∞ so that all the moments R R x k w ( x ) dx exist. Then the Riemann-Hilbert problem(RHP) consists of the following: find a 2 × Y n ( z ) = Y ( z )such that(i) Y ( z ) is analytic for z ∈ C \ R . ISCRETE INTEGRABLE SYSTEMS GENERATED BY HP APPROXIMANTS 9 (ii) Y possesses continuous boundary values for x ∈ R denoted by Y + ( x ) and Y − ( x ), where Y + ( x ) and Y − ( x ) are the limiting values of Y ( z ′ ) as z ′ ap-proaches x from above and below, respectively, and(3.4) Y + ( x ) = Y − ( x ) (cid:18) w ( x )0 1 (cid:19) , x ∈ R . (iii) Y ( z ) has the following asymptotic behavior at infinity:(3.5) Y ( z ) = (cid:18) I + O (cid:18) z (cid:19)(cid:19) (cid:18) z n z − n (cid:19) , z → ∞ . Before giving the solution of this RHP for Y , let us recall that the monic orthogonalpolynomials π n ( z ) = z n + . . . satisfy the following three term recurrence relation: zπ j ( z ) = π j +1 ( z ) + a j π j ( z ) + b j − π j − ( z ) , j ∈ Z + . According to [13], the matrix valued function Y ( z ) given by(3.6) Y ( z ) = π n ( z ) πi R R π n ( x ) w ( x ) x − z dx − πiγ n − π n − ( z ) − γ n − R R π n − ( x ) w ( x ) x − z dx is the unique solution of the RHP for Y . Here γ n is the leading coefficient of thecorresponding orthonormal polynomial. Observe that det Y is an analytic functionon C \ R which has no jump on the real axis. Therefore, det Y is an entire function.Its behaviour near infinity is det Y ( z ) = 1 + O (1 /z ). Thus by Liouville’s theoremwe find that det Y = 1. Consequently, one can consider the matrix W n = Y n +1 Y − n . Clearly W n is an analytic function on C \ R , and since Y n and Y n +1 have the samejump matrix on R we see that L n has no jump on R . Hence W n is an entire matrixfunction. We write the asymptotic condition in the following form Y n ( z ) = (cid:18) I + A ( n ) z + O (1 /z ) (cid:19) (cid:18) z n z − n (cid:19) , where A ( n, m ) is the 2 × /z in the O (1 /z ) term. After somecalculations and using Liouville’s theorem, we find that(3.7) W n = (cid:18) z + A , ( n + 1) − A , ( n ) − A , ( n ) A , ( n + 1) 0 (cid:19) , n ∈ Z + . Remark . In fact we haven’t fully used the Riemann-Hilbert problem to recoverthe wave function and transition matrices in this case. What we actually exploitedis the fact that the solution admits the following factorization Y n ( z ) = R n ( z ) (cid:18) R R w ( x ) dxx − z (cid:19) , where R n ( z ) is a matrix polynomial that has the form R n = W n − . . . W . Basically, R n has a structure similar to that of W [ n, (see formula (3.1)). Moreover, R n coincides with W [ n, up to a constant factor and the inversion. Now, we canclearly see that what we really need here is the Cauchy transform R R w ( x ) dxx − z and itsasymptotic behavior at infinity in order to have (3.5). Therefore, it is clear that one can repeat all the steps for any Borel measure with finite moments of all orders.In other words, we have arrived at the Schur-Euclid algorithm:(i) We start with the function Y ( z ) = Y ( z ) = (cid:18) R R dµ ( x ) x − z (cid:19) , where µ is a probability Borel measure with finite moments of all orders;(ii) Having constructed Y n , we look for the transition matrix L n, of the form(3.7) such that the function Y n +1 = W n Y n obeys the asymptotic condition (3.5).Let us emphasize that the transition matrix in step (ii) is uniquely determined dueto the construction.In the next section we will see that Riemann-Hilbert problems admit gener-alizations in higher dimensions. Thus, they can serve as a tool to develop themultidimensional Schur-Euclid algorithm.4. Hermite-Pad´e and Multiple orthogonal polynomials
Here we present a discrete integrable system associated with a family of Her-mite-Pad´e approximants and multiple orthogonal polynomials.4.1.
Two-dimensional recurrence relations.
We begin by recalling a general-ization of orthogonal polynomials to Hermite-Pad´e polynomials P n,m for two func-tions ( f , f ), which are analytic in a neighbourhood of infinity. It follows fromthe Cauchy theorem applied to (1.6) that Hermite-Pad´e polynomials satisfy theorthogonality relations(4.1) I Γ P n ,n ( z ) z k f j ( z ) dz = 0 , k = 0 , , . . . , n j − , j = 1 , , where the contour Γ := ∂ Ω is the boundary of a domain Ω ∋ ∞ in which thefunctions f j ∈ H (Ω) , j = 1 , f , f ) are the Cauchy transforms(1.12) of positive measures µ , µ with compact support on the real line. In thiscase, the coefficients of the Laurent series (1.5) for ( f , f ) can be considered as themoments of µ , µ : s ( j ) k = I Γ z k f j ( z ) dz −→ s ( j ) k = Z x k dµ j ( x ) , j = 1 , . Using the determinant of the coefficients s ( j ) k (4.2) S n,m = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s (1)0 s (1)1 · · · s (1) n − s (1)1 s (1)2 · · · s (1) n ... ... · · · ... s (1) n + m − s (1) n + m · · · s (1)2 n + m − s (2)0 s (2)1 · · · s (2) m − s (2)1 s (2)2 · · · s (2) m ... ... · · · ... s (2) n + m − s (2) n + m · · · s (2) n +2 m − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , ISCRETE INTEGRABLE SYSTEMS GENERATED BY HP APPROXIMANTS 11 we can write a formula for the Hermite-Pad´e polynomials (4.3) P n,m ( x ) = 1 S n,m (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s (1)0 s (1)1 · · · s (1) n − s (1)1 s (1)2 · · · s (1) n ... ... · · · ... s (1) n + m s (1) n + m +1 · · · s (1)2 n + m − s (2)0 s (2)1 · · · s (2) m − s (2)1 s (2)2 · · · s (2) m ... ... · · · ... s (2) n + m s (2) n + m +1 · · · s (2) n +2 m − x ... x n + m (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) provided that S n,m is nonvanishing. The latter case is a criterion of normality ofthe index ( n, m ). In this paper we assume that all multi-indices are normal and weinvestigate the nearest-neighbor recurrence relations.In [35] a matrix Riemann-Hilbert problem formulation for multiple orthogonalwas proposed. Here we slightly generalize this approach for the case of Hermite-Pad´e polynomials. We can formulate the following Riemann-Hilbert problem: finda 3 × Y such that(i) Y is analytic on C \ Γ , i.e., Y ∈ H (Ω) and Y ∈ H ( C \ Ω),(ii) the continuous limits Y + ( x ) := lim Ω ∋ ξ → x ∈ Γ Y ( ξ ), Y − ( x ) := lim { C \ Ω }∋ ξ → x ∈ Γ Y ( ξ )exist and(4.4) Y + ( x ) = Y − ( x ) f ( x ) f ( x )0 1 00 0 1 , x ∈ Γ , (iii) for z → ∞ one has(4.5) Y ( z ) = (cid:0) I + O (1 /z ) (cid:1) z n + m z − n
00 0 z − m . Following [13], [35] it is easy to show that this Riemann-Hilbert problem has aunique solution in terms of the Hermite-Pad´e polynomials when ( n, m ), ( n − , m )and ( n, m −
1) are normal indices, i.e.,(4.6) Y = P n,m C ( P n,m ) C ( P n,m ) c ( n, m ) P n − ,m c C ( P n − ,m ) c C ( P n − ,m ) c ( n, m ) P n,m − c C ( P n,m − ) c C ( P n,m − ) where the Cauchy transform is used C j ( P ) = 12 πi I Γ P ( x ) f j ( x ) x − z dx, j = 1 , , and the constants c and c are given by − πic ( n, m ) = I Γ P n − ,m ( x ) x n − f ( x ) dx, − πic ( n, m ) = I Γ P n,m − ( x ) x m − f ( x ) dx. One of the natural outcomes of representing the Hermite-Pad´e polynomials in theform of Riemann-Hilbert problems is the nearest-neighbour recurrence relations.
Proposition 4.1.
Suppose all multi-indices ( n, m ) ∈ Z are normal. Then theHermite-Pad´e polynomials satisfy the system of recurrence relations (1.7) .Proof. As we already mentioned in the introduction, the recurrence relations (1.7)for multiple orthogonal polynomials were obtained in [34]. Here, we will followthe same approach. Actually, the proof presented in [34] uses the Riemann-Hilbertproblem but, as we see later, the main ingredient of that proof is a factorization similar to the one revealed in Remark 3.1. Basically, the proof goes along the samelines as the construction of the wave function from the Riemann-Hilbert problemin the case of orthogonal polynomials (see Section 3.2).Let us start by making the standard observation that det Y is an analytic functionin C \ R with no jump on the contour Γ. Hence det Y is an entire function and itsbehavior near infinity is det Y ( z ) = 1 + O (1 /z ). Thus, by Liouville’s theorem wefind that det Y = 1. We can therefore consider the matrix L n,m = Y n +1 ,m Y − n,m , where the subscript ( n, m ) is used for the solution (4.6) of the Riemann-Hilbertproblem with the polynomial P n,m in the entry of the first row and the first columnof Y n,m . Clearly L n,m is an analytic function on C \ R , and since Y n,m and Y n +1 ,m have the same jump matrix on R we see that L n,m has no jump on R . Hence R isan entire matrix function. If we write the asymptotic condition (4.5) as Y n,m ( z ) = (cid:18) I + A ( n, m ) z + O (1 /z ) (cid:19) z n + m z − n
00 0 z − m , where A ( n, m ) is the 3 × /z in the O (1 /z ) term of (4.5),then after some calculus and in view of Liouville’s theorem we find(4.7) L n,m = z + A , ( n + 1 , m ) − A , ( n, m ) − A , ( n, m ) − A , ( n, m ) A , ( n + 1 , m ) 0 0 A , ( n + 1 , m ) 0 1 , where A i,j ( n, m ) is the entry on row i and column j of A ( n, m ). We can thereforewrite(4.8) Y n +1 ,m = L n,m Y n,m . In a similar way we also have(4.9) Y n,m +1 = M n,m Y n,m , with(4.10) M n,m = z + A , ( n, m + 1) − A , ( n, m ) − A , ( n, m ) − A , ( n, m ) A , ( n, m + 1) 1 0 A , ( n, m + 1) 0 0 . Now, introducing(4.11) c n,m = A , ( n, m ) − A , ( n + 1 , m ) , d n,m = A , ( n, m ) − A , ( n, m + 1)and(4.12) a n,m = c ( n, m ) A , ( n, m ) , b n,m = c ( n, m ) A , ( n, m )we see that the (1 , P n +1 ,m ( x ) = ( x − c n,m ) P n,m ( x ) − a n,m P n − ,m ( x ) − b n,m P n,m − ( x ) , and (4.9) gives the second relation in (1.7) P n,m +1 ( x ) = ( x − d n,m ) P n,m ( x ) − a n,m P n − ,m ( x ) − b n,m P n,m − ( x ) . (cid:3) ISCRETE INTEGRABLE SYSTEMS GENERATED BY HP APPROXIMANTS 13
Discrete integrable systems represented by × matrices. Another im-mediate consequence of the reformulation of Hermite-Pad´e approximation in termsof Riemann-Hilbert problems is a bridge between the corresponding recurrence re-lations and discrete integrable system whose transition matrices are 3 × Proposition 4.2.
Let ( f , f ) be a perfect system., i.e., all the multi-indices ( n, m ) are normal. Then there exists a wave function (4.6) on Z and its transitionmatrices are given by (1.8) : L n,m = x + α (1) n,m α (2) n,m α (3) n,m α (4) n +1 ,m α (5) n +1 ,m , M n,m = x + β (1) n,m α (2) n,m α (3) n,m α (4) n,m +1 α (5) n,m +1 , whose entries are related to the coefficients of the recurrence relations (1.7) for theHermite-Pad´e polynomials of the functions f and f as follows: (4.13) c n,m = − α (1) n,m , d n,m = − β (1) n,m , a n,m = − α (4) n,m α (2) n,m , b n,m = − α (5) n,m α (3) n,m . Proof.
We take a function Y of the form (4.6), then (4.8) and (4.9) give us transitionmatrices L n,m , M n,m of the form (4.7) and (4.10). Taking into account that thenormalizing factors in (4.12) are c ( n, m ) = A , ( n, m ) and c ( n, m ) = A , ( n, m ) , the relations (4.12) and (4.11) give (4.13). Finally, we notice that to prove (4.8)and (4.9) we only used the fact that Y admits the following factorization Y ( z ) = R ( z ) H Γ f ( x ) x − z dx H Γ f ( x ) x − z dx , where R is a matrix polynomial. (cid:3) Remark . As we saw in Section 3.1, a continued fraction is just a sequence of2 × × Y ( z ) = Y ( z ) = H Γ f ( x ) x − z dx H Γ f ( x ) x − z dx , where f and f are Laurent series (1.5);(ii) Having constructed Y n,m , we look for the transition matrices L n,m and M n,m of the form (4.7) and (4.10), such that the functions Y n +1 ,m = L n,m Y n,m , Y n,m +1 = M n,m Y n,m , obey the corresponding asymptotic condition (4.5). Note that the transition matrices in step (ii) are uniquely determined due to theconstruction. In Section 5 this idea will be further developed and conventionalcontinued fractions will appear there.Now we simplify the zero curvature condition0 = L n,m +1 M n,m − M n +1 ,m L n,m , to the form of (1.10). Proof of Theorem 1.1 ′ . In [34] the consistency condition for the recurrence coeffi-cients of (1.7) was obtained in the following form:(4.14) d n +1 ,m − d n,m = c n,m +1 − c n,m ,b n +1 ,m − b n,m +1 + a n +1 ,m − a n,m +1 = det d n +1 ,m d n,m c n,m +1 c n,m ! ,a n,m +1 a n,m = c n,m − d n,m c n − ,m − d n − ,m ,b n +1 ,m b n,m = c n,m − d n,m c n,m − − d n,m − . Using the first equation in (4.14), we subtract the columns of the determinantof the second equation in (4.14). We thus obtain the first two equations of (1.10).The third and fourth equations of (1.10) and (4.14) are the same. (cid:3)
Remark . There are other systems related to the concept of orthogonality forwhich the consistency leads to non-trivial zero curvature conditions [20], [29], [30].It turns out that the consistency conditions (4.14) (or, equivalently, the zerocurvature condition) are also sufficient for a sequence of Hermite-Pad´e polynomialsto exist and correspond to a perfect system of functions. To complete the proof ofTheorem 1.1 ′ it remains to prove the following result. Proposition 4.5.
Suppose that the zero curvature condition L n,m +1 M n,m − M n +1 ,m L n,m , holds for a family of invertible matrices L n,m and M n,m of the form (4.7) and (4.10) . Then there are two functions f and f such that the polynomials P n,m satisfying the corresponding relations (1.7) are the Hermite-Pad´e polynomials for f and f .Proof. To determine the functions we first construct the polynomials P n, and P ,m .This can be done since they satisfy ordinary three-term recurrence relations. Sothese polynomials are orthogonal polynomials due to the spectral theorem for or-thogonal polynomials [18, § f be the function corresponding to P n, andlet f be the function for P ,m . Next, the consistency conditions (4.14) allow usto define Y n,m in a unique way for all pairs ( n, m ) ∈ Z . Due to the asymptoticcondition (4.5), the first column of Y n,m consists of Hermite-Pad´e polynomials. Atthe same time, these polynomials coincide with P n,m . Some more details on how toreconstruct the sequence P n,m from the marginal orthogonal polynomials are givenin [12]. (cid:3) ISCRETE INTEGRABLE SYSTEMS GENERATED BY HP APPROXIMANTS 15
To conclude this subsection, note that the wave function Ψ n,m coincides with Y n,m for ( n, m ) ∈ Z and it can be extended to the entire lattice Z by the symmetry(4.15) Ψ − n,m = Ψ n,m , Ψ n, − m = Ψ n,m , Ψ − n, − m = Ψ n,m , n, m ∈ Z + . The underlying boundary value problem
In this section we discuss the discrete system and give a few algorithms of recon-structing solutions from boundary data. Finally we will consider some observableclasses of initial data for which the system is solvable.Here we are concerned with the question of finding a solution of the differenceequations(5.1) c n,m +1 = c n,m + ( a + b ) n +1 ,m − ( a + b ) n,m +1 ( c − d ) n,m ,d n +1 ,m = d n,m + ( a + b ) n +1 ,m − ( a + b ) n,m +1 ( c − d ) n,m ,a n,m +1 = a n,m ( c − d ) n,m ( c − d ) n − ,m ,b n +1 ,m = b n,m ( c − d ) n,m ( c − d ) n,m − , n, m ≥ , subject to the boundary conditions a ,m = 0 , d ,m = ˆ d m , b ,m +1 = ˆ b m +1 , m ∈ Z + ,b n, = 0 , c n, = ˆ c n , a n +1 , = ˆ a n +1 , n ∈ Z + , (5.2)where ˆ c n , ˆ d m , ˆ a n +1 , and ˆ b m +1 are sequences of complex numbers. More precisely,we are interested in solutions for which one has a n,m = 0 , b n,m = 0 , n, m > . According to Theorem 1.1 ′ , such a solution exists if and only if there is a perfectsystem of two functions f and f . In this case the initial data ˆ c n , ˆ a n +1 , ˆ d m , andˆ b m +1 are the entries of the J -fraction representations of f and f : f ( z ) ∼ − z − ˆ c − ˆ a z − ˆ c − ˆ a . . . , f ( z ) ∼ − z − ˆ d − ˆ b z − ˆ d − ˆ b . . . , where we have(5.3) a n +1 , = ˆ a n +1 = 0 , b ,m +1 = ˆ b m +1 = 0 . Basically, this means that in order to have a well-posed boundary value problem,the initial data have to be generated through the Schur-Euclid algorithm by twofunctions that form a perfect system, i.e., a system of two functions that determinesthe entire table of multiple orthogonal polynomials.One of the main obstacles to construct a table of multiple orthogonal polynomialsis to ensure that each index is normal, that is, the corresponding determinant S n,m is non-zero. It cannot be done for two arbitrary analytic functions or evenfor any two positive measures and this issue was addressed for the first timeby K. Mahler [21], who coined the notion of perfect systems. To be more precise,a system of two measures is called perfect if each index in the corresponding table is normal, i.e., S n,m = 0 for all n, m ∈ Z + . At the end of this section we givetwo rather general classes of perfect systems. In turn, these systems give rise to aninfinite number of solutions of the discrete integrable system in question. Beforegoing into more details about those classes, we reformulate a part of Theorem 1.1 ′ and give a constructive proof, which actually is the way to solve (5.1), (5.2). Proposition 5.1.
If the initial data ˆ c n , ˆ d m , ˆ a n +1 , and ˆ b m +1 correspond to a perfectsystem, then the boundary value problem (5.1) , (5.2) has a solution that satisfiesthe condition (5.3) . Moreover, the problem can be solved in the following way.The boundary data ˆ c n , ˆ d m , ˆ a n +1 , and ˆ b m +1 define the moments and the solution ( c, a ) n,m , ( d, b ) n,m can be recovered from the moments.Proof. The statement is a straightforward consequence of Theorem 1.1 ′ . So, weknow that the underlying system of functions f and f is perfect, that is, we havethat the corresponding moments are such that S n,m = 0 for all n, m ∈ Z + . Asa matter of fact, it is a standard technique that recovers the moments from theentries of the corresponding J -fraction (see [2]). Furthermore, using the momentsone can reconstruct the solution in the following manner. The coefficients a n,m and b n,m can be found via the formulas from [34] a n,m = S n +1 ,m S n − ,m (cid:0) S n,m (cid:1) , b n,m = S n,m +1 S n,m − (cid:0) S n,m (cid:1) . We also know from [34] that(5.4) d n,m − c n,m = S n,m S n +1 ,m +1 S n +1 ,m S n,m +1 . Finally, the rest are found by summation of the first and second equations of thesystem for consecutive indices c n,m +1 = c n, + m X i =1 ( a + b ) n +1 ,i − ( a + b ) n,i +1 ( c − d ) n,i ,d n +1 ,m = d ,m + n X i =1 ( a + b ) i +1 ,m − ( a + b ) i,m +1 ( c − d ) i,m , since the elements on the right hand sides are already determined or are part of theinitial data. (cid:3) Sometimes, given a perfect system, it is easier to find the Hermite-Pad´e poly-nomials rather then moments and for this reason we can elaborate a bit more onthe idea mentioned in Remark 4.3 in order to see some continued fractions in theconventional sense. Following [14] we introduce the following functions m (1) − ( z, n, m ) = − P n,m ( z ) P n +1 ,m ( z ) , m (2) − ( z, n, m ) = − P n,m ( z ) P n,m +1 ( z ) , which can serve as a tool for solving the system (1.10). Proposition 5.2.
If the initial data ˆ c n , ˆ d m , ˆ a n +1 , and ˆ b m +1 correspond to a perfectsystem, then the boundary value problem (5.1) , (5.2) has a solution that satisfies thecondition (5.3) . Furthermore, the boundary value problem (5.1) , (5.2) can be solvedin the following way. From the initial data ( c, a ) n, , ( d, b ) ,m one can reconstructthe family of polynomials P n,m (which can formally be done via formula (4.3) ). In ISCRETE INTEGRABLE SYSTEMS GENERATED BY HP APPROXIMANTS 17 turn, these polynomials define the functions m (1) − ( z, n, m ) and m (2) − ( z, n, m ) , whichadmit the following continued fraction expansions m (1) − ( z, n, m ) = − z − c n,m − a n,m z − c n − ,m − . . . − b n,m z − d n,m − − . . . ,m (2) − ( z, n, m ) = − z − d n,m − a n,m z − c n − ,m − . . . − b n,m z − d n,m − − . . . that determine the solution ( c, a ) n,m , ( d, b ) n,m to the equation (1.10) .Proof. It is not so hard to see that the relations (1.7) can be rewritten as thefollowing generalization of the discrete Riccati equation: m (1) − ( z, n, m ) = − z − c n,m + a n,m m (1) − ( z, n − , m ) + b n,m m (2) − ( z, n, m − , (5.5) m (2) − ( z, n, m ) = − z − d n,m + a n,m m (1) − ( z, n − , m ) + b n,m m (2) − ( z, n, m − . (5.6)Then, we see that applying (5.5) and (5.6) iteratively leads to the continued frac-tion expansions under consideration, which allow us to recursively reconstruct thesolution from the initial boundary data. Namely, it is clear how to reconstruct c n,m and d n,m for all indices from m (1) − ( z, n, m ) and m (2) − ( z, n, m ) in the first place. Thenthe first term of the asymptotic expression − m (1) − ( z, n, m ) − z + c n,m , z → ∞ (or, equivalently, the analogous one for m (2) − ( z, n, m )) determines a n,m + b n,m = f n,m and the consecutive term gives a n,m c n − ,m + b n,m d n,m − = g n,m . Next, taking intoaccount the first equation in (4.14) we have that c n − ,m − d n,m − = c n − ,m − − d n − ,m − . Since the system is perfect we get that c n − ,m − − d n − ,m − = 0 (see (5.4)). There-fore, a n,m and b n,m are uniquely determined. (cid:3) Let us emphasize here that the above given continued fractions branch into twocontinued fractions on each level. Actually, what we see is that there are twofractions behind the scene. On the one hand, they are similar to the classical onein (3.3) but, on the other hand, they have a certain two-dimensional structure andthe two fractions are identical except for one entry.
Remark . There is one more algorithm available, which is obtained by combiningthe well-known Jacobi-Perron algorithm (see [24] for the details) and a result from[12]. Namely, the Jacobi-Perron algorithm expands the vector ( f , f ), where f and f form a perfect system, into a vector continued fraction. The approximantsof this vector continued fraction consists of rational functions with the same de-nominator. More importantly, the denominators are the Hermite-Pad´e polynomialsthat correspond to the so-called step-line, that is when m = n and m = n −
1. Fur-thermore, the step-line polynomials satisfy recurrence relations whose coefficientsare the input for the algorithm given in [12] to reconstruct all the coefficients of thenearest neighbour recurrence relations from the step-line.
Angelesco systems.
A. Angelesco considered in [3] the following systems ofmeasures. Let ∆ and ∆ be disjoint bounded intervals on the real line and ( µ , µ )be a system of measures such that supp µ j = ∆ j .Fix ~n = ( n , n ) ∈ Z and consider the multiple orthogonal polynomials ofthe so called Angelesco system ( b µ , b µ ) relative to ~n . Here b µ denotes the Cauchytransform of µ : b µ ( z ) = Z dµ ( x ) z − x . By construction, we have that Z x ν P ~n ( x ) dµ j ( x ) = 0 , ν = 0 , . . . , n j − , j = 1 , . Therefore, P ~n has n j simple zeros in the interior (with respect to the Euclideantopology of R ) of ∆ j . As a consequence, since the intervals ∆ j are disjoint, deg P ~n = | ~n | = n + n and Angelesco systems are perfect .5.2.
Nikishin systems.
Unfortunately, Angelesco’s paper received little attentionand such systems reappeared many years later in [22] where E.M. Nikishin deducedsome of their formal properties. He also introduced another class of systems forwhich the perfectness was proved only recently [23].To get an idea about these systems, let us consider two disjoint bounded intervals∆ , ∆ on the real line. Suppose we are given two measures σ and σ supportedon ∆ and ∆ , respectively. With these two measures we define a third one in thefollowing way d h σ , σ i ( x ) = b σ ( x ) dσ ( x );that is, one multiplies the first measure by a weight formed by the Cauchy transformof the second measure. Thus, we have arrived at the notion of a Nikishin system.A system of two measures ( µ , µ ) of the form dµ ( x ) = dσ ( x ) , dµ ( x ) = Z dσ ( t ) x − t dµ ( x ) = d h σ , σ i is called a Nikishin system (of order 2). In fact, a Nikishin system can be definedfor any finite number of measures. We also emphasize that the measures from aNikishin system have the same support, which is a totally different situation thanin the case of an Angelesco system.Finally, it is worth mentioning here that it was a long standing problem to provethat a general Nikishin system ( µ , µ , . . . , µ p ) is perfect for p ≥
2. This fact wasfinally proved in the remarkable paper [11].
References [1] V. E. Adler,
Discrete equations on planar graphs. Symmetries and integrability of differenceequations (Tokyo, 2000) , J. Phys. A (2001), no. 48, 10453–10460.[2] N. I. Akhiezer, The Classical Moment Problem and Some Related Problems in Analysis ,Hafner publishing Co., New York, 1965.[3] A. Angelesco,
Sur deux extensions des fractions continues alg´ebriques , C.R. Acad. Sci. Paris (1919), 262–263.[4] R. Ap´ery, Irrationalit´e de ζ (2) et ζ (3), Ast´erisque (1979), 11–13.[5] A. I. Aptekarev, Multiple orthogonal polynomials , J. Comput. Appl. Math. (1998), no. 1–1,423–448.[6] A. I. Aptekarev, Rational approximants for Euler’s constant and recurrence relations , Pro-ceedings of the Steklov Institute of Mathematics (2011), 138–141.
ISCRETE INTEGRABLE SYSTEMS GENERATED BY HP APPROXIMANTS 19 [7] A. I. Aptekarev, A. B. J. Kuijlaars,
Hermite-Pad´e approximations and multiple orthogonalpolynomial ensembles , Russ. Math. Surv. (2011), no. 6, 1133–1200.[8] A. I. Bobenko, Discrete differential geometry. Integrability as consistency , Discrete integrablesystems, 85–110, Lecture Notes in Phys. , Springer, Berlin, 2004.[9] A. I. Bobenko, Yu. B. Suris,
Integrable systems on quad-graphs , Int. Math. Res. Not. 2002,no. 11, 573–611.[10] A. I. Bobenko, Yu. B. Suris,
Discrete Differential Geometry: Integrable Structure , GraduateStudies in Mathematics, Vol. 98, AMS, 2008.[11] U. Fidalgo Prieto, G. L´opez Lagomasino,
Nikishin systems are perfect , Constr. Approx. (2011), no. 3, 297–356.[12] G. Filipuk, M. Haneczok, W. Van Assche, Computing recurrence coefficients of multipleorthogonal polynomials , arXiv:1406.0364.[13] A.S. Fokas, A.R. Its, A.V. Kitaev,
The isomonodromy approach to matrix models in 2Dquantum gravity , Comm. Math. Phys. (1992), no. 2, 395–430.[14] F. Gesztesy, B. Simon, m-functions and inverse spectral analysis for finite and semi-infiniteJacobi matrices . J. Anal. Math. 73 (1997), 267–297.[15] C. Hermite,
Sur la fonction exponentielle , C.R. Acad. Sci. Paris (1873), 18–24; 74–79;226–233.[16] A. Ya. Khinchin, Continued fractions . Translated from the third (1961) Russian edition.Reprint of the 1964 translation. Dover Publications, Inc., Mineola, NY, 1997.[17] J. van Iseghem,
Vector Stieltjes continued fraction and vector QD algorithm . InternationalConference on Numerical Algorithms, Vol. I (Marrakesh, 2001). Numer. Algorithms 33 (2003),no. 1-4, 485–498.[18] M.E.H. Ismail,
Classical and Quantum Orthogonal Polynomials in One Variable , Encyclo-pedia of Mathematics and its Applications , Cambridge University Press, 2005.[19] A. B. J. Kuijlaars, Multiple orthogonal polynomials in random matrix theory , Proceedingsof the International Congress of Mathematicians. Volume III, 1417–1432, Hindustan BookAgency, New Delhi, 2010.[20] V. Papageorgiou, B. Grammaticos, A. Ramani,
Orthogonal polynomial approach to discreteLax pairs for initial-boundary value problems of the QD algorithm , Lett. Math. Phys. (1995), no. 2, 91–101.[21] K. Mahler, Perfect systems , Compos. Math. (1968), no. 2, 95–166.[22] E. M. Nikishin, A system of Markov functions , Vestnik Moskov. Univ. Ser. I Mat. Mekh(1979):4, 60–63 (in Russian); translation in Moscow Univ. Math. Bull. (1979), 63–66.[23] E. M. Nikishin, On simultaneous Pad´e approximants , Matem. Sb. (1980), 499–519 (inRussian); translation in Math. USSR Sb. (1982), 409–425.[24] E. M. Nikishin, V. N. Sorokin, Rational approximations and orthogonality . Translated fromthe Russian by Ralph P. Boas. Translations of Mathematical Monographs, 92. AmericanMathematical Society, Providence, RI, 1991.[25] S. P. Novikov,
Four Lectures: Discretization and Integrability. Discrete Spectral Symmetries ,in “Integrability”, Lecture Notes in Physics , Springer, Berlin, 2009, pp. 119–138.[26] L. A. Sakhnovich,
Interpolation theory and its applications , Mathematics and its Applications , Kluwer Academic Publishers, Dordrecht, 1997.[27] L. A. Sakhnovich,
Spectral theory of canonical differential systems. Method of operator iden-tities , Operator Theory: Advances and Applications , Birkh¨auser Verlag, Basel, 1999.[28] S. Smirnov,
Discrete complex analysis and probability , Proceedings of the International Con-gress of Mathematicians. Volume I, 595–621, Hindustan Book Agency, New Delhi, 2010.[29] P. E. Spicer, F. W. Nijhoff, P. H. van der Kamp,
Higher analogues of the discrete-time Todaequation and the quotient-difference algorithm , Nonlinearity (2011), no. 8, 2229–2263.[30] V. Spiridonov, A. Zhedanov, Spectral transformation chains and some new biorthogonal ra-tional functions , Comm. Math. Phys. (2000), no. 1, 49–83.[31] Y. B. Suris,
Bi-Hamiltonian structure of the qd algorithm and new discretizations of theToda lattice , Phys. Lett. A (1995), no. 3–4, 153–161.[32] W. Van Assche, Multiple orthogonal polynomials, irrationality and transcendence , in “Con-tinued fractions: from analytic number theory to constructive approximation”, ContemporaryMathematics (1999), 325-342.[33] W. Van Assche,
Little q -Legendre polynomials and irrationality of certain Lambert series ,Ramanujan J. (2001), 295–310. [34] W. Van Assche, Nearest neighbor recurrence relations for multiple orthogonal polynomials ,J. Approx. Theory (2011), no. 10, 1427–1448.[35] W. Van Assche, J. S. Geronimo, A. B. J. Kuijlaars,
Riemann-Hilbert problems for multipleorthogonal polynomials , Special functions 2000: current perspective and future directions(Tempe, AZ), NATO Sci. Ser. II Math. Phys. Chem. , Kluwer Acad. Publ., Dordrecht,2001, pp. 23–59. Alexander I. Aptekarev, Keldysh Institute for Applied Mathematics, Russian Acad-emy of Sciences, Miusskaya pl. 4, 125047 Moscow, RUSSIAMaxim Derevyagin, University of Mississippi, Department of Mathematics, Hume Hall305, P. O. Box 1848, University, MS 38677-1848, USA
E-mail address : [email protected]@gmail.com