Multiple orthogonal polynomials, $d$-orthogonal polynomials, production matrices, and branched continued fractions

I analyze an unexpected connection between multiple orthogonal polynomials, $d$-orthogonal polynomials, production matrices and branched continued fractions. This work can be viewed as a partial extension of Viennot's combinatorial theory of orthogonal polynomials to the case where the production matrix is lower-Hessenberg but is not necessarily tridiagonal.


Introduction
The goal of this paper is to point out, and then analyze in detail, an unexpected connection between multiple orthogonal polynomials and d-orthogonal polynomials on the one hand, and production matrices and branched continued fractions on the other -objects that arose over the past few decades in the special-functions and enumerative-combinatorics communities, respectively.It is appropriate to begin, therefore, by explaining briefly each of these four concepts.
Closely related to multiple orthogonal polynomials are the so-called d-orthogonal polynomials.A sequence (P n (x)) n≥0 of monic polynomials is said to be d-orthogonal [50,77] with respect to a sequence Γ 0 , . . ., Γ d−1 of linear forms in case Γ k (x P n (x)) = 0 whenever n > d +k.(To avoid trivialities, it is also usually required that Γ k (x P n (x)) = 0 when n = d + k.)For d = 1 this reduces to the ordinary concept of orthogonality.It turns out that the sequence of multiple orthogonal polynomials of type II taken along the so-called stepline is d-orthogonal (for d = r) with respect to the linear forms associated to the measures µ 1 , . . ., µ r .
Production matrices [23,24] have become in recent years an important tool in enumerative combinatorics (see Section 2.2 for a brief summary).In the special case of a tridiagonal production matrix, this construction goes back to Stieltjes' [70,71] work on continued fractions: the production matrix of a classical S-fraction or J-fraction is tridiagonal.Moreover, the classical J-fraction and tridiagonal production matrix associated to the moment sequence of a measure µ are closely related to the sequence of orthogonal polynomials associated to µ.This connection was comprehensively investigated by Viennot [80,81] in the early 1980s, who developed a general combinatorial theory of orthogonal polynomials, building on Flajolet's [31] combinatorial theory of continued fractions.Our work here can be viewed as a partial extension of Viennot's theory to the case where the production matrix is lower-Hessenberg (i.e.vanishes above the first superdiagonal) but is not necessarily tridiagonal.Indeed, this extension was already begun by Viennot himself [80, sections III.5 and V.6].
Various types of branched continued fractions have been introduced in the analysis literature [10,11,22] [ 49, pp. 274-280, 285] [20, p. 28], but we are not concerned here with these.Rather, we are concerned with the branched continued fractions that have been introduced by combinatorialists and whose Taylor coefficients are the generating polynomials for selected types of lattice paths, generalizing the work of Flajolet [31] on classical continued fractions.This investigation was also initiated by Viennot [80, section V.6], who briefly considered the branched continued fractions (fractions multicontinuées) generated by Lukasiewicz paths.This work was carried forward in the Ph.D. theses of Roblet [59] and Varvak [78].Further applications were made by Gouyou-Beauchamps [33] and Drake [28].Subsequently, Albenque and Bouttier [4] introduced the branched continued fractions generated by m-Dyck paths and proved many interesting results about them.Most recently, Pétréolle, Sokal and Zhu [56] carried out a comprehensive analysis of the branched continued fractions associated to m-Dyck, m-Schröder and Lukasiewicz paths, with emphasis on questions related to total positivity; see also [21,55] for further applications.For these branched continued fractions, the production matrix is lower-Hessenberg but not (except in the classical cases) tridiagonal.See Section 2.4 for a brief summary.
Finally, let us mention the remarkable Ph.D. thesis of Drake [27], who initiated the combinatorial theory of multiple orthogonal polynomials and who foresaw the link with branched continued fractions [27, p. 1].Our work here can be viewed as an extension, and to some extent a completion, of his.
Throughout this paper, we fix a commutative ring (with identity element 1 = 0) R: we will use sequences and matrices with entries in R, and polynomials and formal power series with coefficients in R. The analyst reader should feel free to imagine, without too much loss of generality, that R = R.However, in applications of this formalism there will often be parameters, and we will usually prefer to treat these parameters as algebraic indeterminates ξ; then R will be either the ring R[ξ] of polynomials in these indeterminates or the field R(ξ) of rational functions in these indeterminates.
In particular, to any (positive or signed) measure µ on R that has finite moments of all orders, there is canonically associated a linear functional L on the polynomial ring R[x], defined by L(x n ) = x n dµ(x).It is well known [14] that the theory of orthogonal polynomials (or at least the simplest parts of it) can be expressed entirely in terms of this linear functional -or equivalently, in terms of the sequence of moments a n = L(x n ) -without reference to the measure µ.We shall adopt this approach here, and also replace R by an arbitrary commutative ring R.
The plan of this paper is as follows: In Section 2 we collect some basic definitions and facts concerning multiple orthogonal polynomials, production matrices and branched continued fractions.In Section 3 we use the theory of production matrices to demonstrate some very simple relations between sequences of monic polynomials, the linear recurrences they satisfy, and their dual sequences of linear functionals.We also generalize Viennot's [80] formula for the expectation of products of orthogonal polynomials.In Section 4 we analyze sequences of monic polynomials that are orthogonal to a sequence of linear functionals.In Section 5 we show how this theory applies to ordinary orthogonal polynomials, and in Section 6 we apply it to multiple orthogonal polynomials.Finally, in Section 7 we examine some concrete examples.In the Appendix we prove a basic result concerning LU factorization for matrices over a commutative ring.

Preliminaries
In this section we provide a brief introduction to multiple orthogonal polynomials [5,51,75] [36,Chapter 23], production matrices [23,24,56,65], and branched continued fractions [56].The reader familiar with one or more of these topics can skim those parts quickly, with the main aim of fixing the notation.

Multiple orthogonal polynomials
We begin by giving a brief introduction to the theory of multiple orthogonal polynomials, following [36,Chapter 23] but making a few comments about "algebraizing" the theory to allow coefficients in an arbitrary commutative ring R. We limit attention to the multiple orthogonal polynomials of type II, since these are the only ones that will arise in the remainder of the paper.(I leave it to others to investigate whether there is any analogue of the connections discussed here for the multiple orthogonal polynomials of type I.) Fix an integer r ≥ 1, and let µ 1 , . . ., µ r be positive measures on the real line with finite moments of all orders.We use multi-indices n = (n 1 , . . ., n r ) ∈ N r and write |n| = n 1 + . . .+ n r .The multiple orthogonal polynomial of type II for the multi-index n is the degree-|n| monic polynomial P n (x) = x |n| + . . .satisfying the orthogonality relations whenever such a polynomial exists and is unique.The equations (2.1) give a system of |n| linear equations for the |n| non-leading coefficients of the polynomial P n (x); the multi-index n is said to be normal whenever the solution exists and is unique.Note that the coefficient matrix of this system is the transpose of the |n| × |n| matrix M,N is the M × N Hankel matrix of the moments of µ j : that is, H M,N = (m Therefore, the multi-index n is normal if and only if det M n = 0. A system of measures µ 1 , . . ., µ r is said to be perfect in case all n ∈ N r are normal.Several general sufficient conditions for a system to be perfect are known (Angelesco systems, AT systems, Nikishin systems, . . .): see [36,Chapter 23] [75].We shall henceforth restrict attention to perfect systems.
The orthogonality conditions (2.1) can be trivially re-expressed in terms of the linear forms L (1) , . . ., L (r) associated to the measures µ 1 , . . ., µ r , which are defined by k by L (j) (x k ).Moreover, from this point of view, the measures µ j need not be positive measures; indeed, the linear forms L (j) need not come from (signed) measures at all.Provided that one can show, one way or another, that all n ∈ N r are normal, the multiple orthogonal polynomials are well-defined.
Having done this, we can go farther and "algebraize" the theory by considering polynomials with coefficients in an arbitrary commutative ring (with identity element 1 = 0) R, rather than just R = R.We fix linear forms L (1) , . . ., L (r) on the polynomial ring R[x], and define "moments" m (j) k = L (j) (x k ); then the multi-index n is normal if and only if det M n is an invertible element of the ring R.
Let us now make a simple but important observation.Fix a multi-index n = (n 1 , . . ., n r ) ∈ N r , and suppose that the polynomial P n (x) satisfies the orthogonality relations (2.1) with respect to some family of (not necessarily positive) measures µ = (µ 1 , . . ., µ r ).Then P n (x) also satisfies the orthogonality relations (2.1) with respect to any family of (not necessarily positive) measures µ = (µ 1 , . . ., µ r ) where µ i is any linear combination of {µ j : n j ≥ n i }.In particular, if n 1 ≥ n 2 ≥ . . .≥ n r , then we can take µ i = i j=1 c ij µ j for any lower-triangular matrix C = (c ij ) 1≤i,j≤r .That is, µ i is an arbitrary linear combination of µ 1 , . . ., µ i .This observation will play an important role in what follows (see Section 4 ff.).
The collection (P n (x)) n∈N r of (monic) multiple orthogonal polynomials of type II satisfies a variety of recurrences, generalizing the well-known three-term recurrence for conventional orthogonal polynomials.Here is one [36,Theorem 23.1.7

et seq.]:
We denote by e k the multi-index with entry 1 in position k and 0 elsewhere.For a permutation π of {1, . . ., r}, we write s (π) j = e π(1) + e π(2) + . . .+ e π(j) for 1 ≤ j ≤ r.Then there exist real numbers a with the convention that P n−s (π) j (x) = 0 whenever one or more of the entries in n−s (π) j is negative.(Note that the coefficients a (π) n,j do not depend on k, and the coefficients a (k) n,0 do not depend on π.But we will never use this fact.)Now let j 1 , j 2 , . . .be an infinite sequence of elements of {1, . . ., r}, and define a sequence (n k ) k≥0 of multi-indices in N r by n k = k i=1 e j i .These multi-indices satisfy |n k | = k and describe an increasing nearest-neighbor path in N r in which the ith step is along direction j i .Now let P k (x) def = P n k (x) be the multiple orthogonal polynomial of type II along this path in N r .It then follows from (2.4) that the singly-indexed sequence ( P n (x)) n≥0 satisfies an (r + 2)-term recurrence of the form where π n,n+1 = 1 and π nk = 0 for k < 0; of course the coefficients π nk depend on the choice of nearest-neighbor path.The recurrence (2.5) will play a central role in the remainder of this paper.Please observe that the coefficients π nk in this recurrence can be collected into a matrix Π = (π nk ) n,k≥0 that is unit-lower-Hessenberg and (r, 1)banded: that is, A particularly important role is played by the multi-indices n = (n 1 , . . ., n r ) lying on the stepline: this is the near-diagonal sequence starting at (0, 0, . . ., 0) and following the path (n, n, . . ., n) → (n + 1, n, . . ., n) → (n + 1, n + 1, . . ., n) → . . .→ (n + 1, n + 1, . . ., n + 1) → . . . .In other words, we define a singly-indexed sequence ( P n (x)) n≥0 by P n (x) = P (n 1 ,...,nr) (x) where The stepline polynomials ( P n (x)) n≥0 are a special case of the nearest-neighbor-path polynomials ( P n (x)) n≥0 , so they satisfy an (r + 2)-term recurrence of the form (2.5).

Production matrices
In this subsection we give a brief introduction to the theory of production matrices [23,24]; see also [65] [56, sections 8.1 and 9.2] for further discussion.In the general theory, the production matrix can be any row-finite or column-finite matrix.
Here, however, we shall restrict attention to production matrices that are unit-lower-Hessenberg.Also, in the general theory the production matrix is usually called P ; but here we shall call it Π in order to avoid confusion with the sequence of polynomials P n (x).
So let Π = (π ij ) i,j≥0 be a unit-lower-Hessenberg matrix (indexed by N) with entries in a commutative ring R: that is, π i,i+1 = 1 and π ij = 0 for j > i + 1.Then let A = (a nk ) n,k≥0 be the matrix defined by a nk = (Π n ) 0k .It is easy to see that A is unit-lower-triangular, i.e. a nn = 1 and a nk = 0 for k > n.Writing out the matrix multiplications explicitly, we have so that a nk is the total weight for all n-step walks in N from i 0 = 0 to i n = k, in which the weight of a walk is the product of the weights of its steps, and a step from i to j gets a weight π ij .(Since Π is lower-Hessenberg, these are Lukasiewicz walks, i.e. the allowed steps are i → j with 0 ≤ j ≤ i + 1.) Yet another equivalent formulation is to define the entries a nk by the recurrence with the initial condition a 0k = δ 0k .We shall call Π the production matrix and A the output matrix , and we write A = O(Π).These definitions can be given a compact matrix formulation.Let ∆ be the matrix with 1 on the superdiagonal and 0 elsewhere, i.e. ∆ n,n+1 = 1 and ∆ nk = 0 for k = n+1 (of course it is unit-lower-Hessenberg).Then for any matrix M with rows indexed by N, the product ∆M is simply M with its zeroth row removed and all other rows shifted upwards.(Some authors use the notation M def = ∆M .)The recurrence (2.8) can then be written as ∆A = AΠ . (2.9) Since A is unit-lower-triangular, it is invertible, so (2.9) is equivalent to Conversely, it is not difficult to see that for any unit-lower-triangular matrix A, the matrix A −1 ∆A is unit-lower-Hessenberg.It therefore follows that for each unit-lowertriangular matrix A, there is a unique unit-lower-Hessenberg matrix Π such that A = O(Π), and it is given by Π = A −1 ∆A.Now let B = A −1 be the inverse of A (which is of course also unit-lower-triangular).We then have a one-to-one correspondence between unit-lower-Hessenberg matrices Π, unit-lower-triangular matrices A and unit-lower-triangular matrices B, defined by In Section 3 we will see how the matrices A, B and Π arise in different characterizations of sequences of monic polynomials.
Remarks. 1. Production matrices are nowadays widely used in enumerative combinatorics: thus, for instance, the entry for a triangular array in the On-Line Encyclopedia of Integer Sequences [53] often gives its production matrix.
2. Several subclasses of lower-Hessenberg production matrices are of especial combinatorial interest: • Tridiagonal production matrices correspond to Motzkin walks (i.e. the allowed steps are i → j with j ∈ {i − 1, i, i + 1}) and thence to classical J-fractions [31], as will be explained in Section 2.3.
• Lower-Hessenberg production matrices of the form 3. When the commutative ring R is equipped with a partial order, production matrices are also a powerful tool for attacking problems related to total positivity [65].In particular, the total positivity of the production matrix Π is a sufficient (but far from necessary) condition for the total positivity of its output matrix O(Π) and for the Hankel-total positivity of the zeroth-column sequence of O(Π).See [65] [56, sections 8.1 and 9.2] [21, 55,64] for precise statements, proofs, and further discussion and applications.

Classical continued fractions
As preparation for the discussion of branched continued fractions in Section 2.4, as well as for some applications later in this paper, we begin by giving a very brief review of selected aspects of the theory of classical continued fractions (J-fractions and S-fractions).We will follow the notation and terminology used nowadays by combinatorialists [31], as this is the most appropriate for our work; but we will also point out the translation to the formalism employed in the classical analysis literature on continued fractions [20,38,49,54,82] and the moment problem [1,3,60,63,71].
We shall consider continued fractions of either Stieltjes (S) type, or Jacobi (J) type, Here these expressions are to be interpreted as formal power series in the indeterminate t; we do not wish to address questions of convergence.Thus, the continuedfraction coefficients α = (α n ) n≥1 , β = (β n ) n≥1 and γ = (γ n ) n≥0 are sequences in a commutative ring R, and they determine the sequence a = (a n ) n≥0 of Taylor coefficients by formal expansion of the continued fraction.Indeed, it is conceptually simplest to consider α, β, γ as algebraic indeterminates; then the a n are polynomials with integer coefficients in these indeterminates: We call S n (α) the Stieltjes-Rogers polynomials, and J n (β, γ) the Jacobi-Rogers polynomials.
In a seminal 1980 paper, Flajolet [31] gave a combinatorial interpretation of the Stieltjes-Rogers and Jacobi-Rogers polynomials in terms of lattice paths.We recall that a Motzkin path of length n is a path in the right quadrant N × N, starting at (0, 0) and ending at (n, 0), using steps (1, 1) ["rise"], (1, 0) ["level step"] and (1, −1) ["fall"].More generally, a Motzkin path at level k is a path in N × N ≥k , starting at (0, k) and ending at (n, k), using the same steps.A Motzkin path is called a Dyck path if it has no level steps; obviously a Dyck path must have even length.
(a) The Jacobi-Rogers polynomial J n (β, γ) is the generating polynomial for Motzkin paths of length n, in which each rise gets weight 1, each level step at height i gets weight γ i , and each fall from height i gets weight β i .
(b) The Stieltjes-Rogers polynomial S n (α) is the generating polynomial for Dyck paths of length 2n, in which each rise gets weight 1 and each fall from height i gets weight α i .
Proof [31].(a) For each k ≥ 0, let f k (t) be the generating function for Motzkin paths at level k (of arbitrary length) in which each rise gets weight 1, each level step at height i gets weight γ i , each fall from height i gets weight β i , and each step of any kind gets an additional weight t.It is a formal power series in the indeterminate t, with coefficients that are polynomials in β, γ.Now let P be any Motzkin path at level k; and if it is of nonzero length, split it at its first return to height k, yielding P = P P .Then P is either a single level step at height k, or else a path of the form U P k+1 D where U is a rise k → k + 1, P k+1 is an arbitrary Motzkin path at level k + 1, and D is a fall k + 1 → k.Furthermore, P is an arbitrary Motzkin path at level k.We thus deduce the functional equation or equivalently . (2.17) Iterating (2.17), we see immediately that f k is given by the continued fraction and in particular that f 0 is given by (2.15).(b) This follows from part (a) by setting γ = 0, renaming β as α, and renaming t 2 as t.
Remarks. 1.In the function-theoretic literature on the moment problem [1,3,60,63,71] and continued fractions [20,38,49,54,82], the generating function for a sequence a = (a n ) n≥0 of real numbers is most often written in the form This formulation has the property that if a is a moment sequence with representing measure µ, i.e.
is analytic in the upper half-plane Im z > 0 and has the series (2.19) as its large-z asymptotic expansion, uniformly in each sector Given a power series of the form (2.19), the S-type continued fraction is then written in the form [63, p. viii] [82, p. 329] (note that l n is multiplied by z for n odd but not for n even), which is easily seen to be equivalent to (2.12) if we normalize to a 0 = 1 (hence l 1 = 1) and then set α 1 = 1/l 2 and α n = 1/(l n l n+1 ) for n ≥ 2; the reverse translation is Likewise, the J-type continued fraction is written in the form [63, pp. viii, 31] which is easily seen to be equivalent to (2.13) if we normalize to λ 1 = 1 and then set γ n = c n+1 and β n = λ n+1 .
2. My use of the terms "S-fraction" and "J-fraction" follows the general practice in the combinatorial literature, starting with Flajolet [31].The classical literature on continued fractions [20,38,49,54,82] generally uses a different terminology.For instance, Jones and Thron [38, pp. 128-129, 386-389] use the term "regular C-fraction" for (a minor variant of) what I have called an S-fraction; they call it an "S-fraction" if all α n < 0. They use the term "associated continued fraction" for (a minor variant of) what I have called a J-fraction, and use the term "J-fraction" for (2.23) with λ n = 0.
3. It is worth observing that an S-fraction can always be transformed into a Jfraction by contraction [82, p. 21] [80, p. V-31]: namely, (2.12) and (2.13) are equal if See [82, pp. 20-22] for the classic algebraic proof; see [30,  Let us now generalize these definitions; we concentrate on the case of J-fractions, but similar constructions can be applied to S-fractions.A partial Motzkin path of length n is a path in the right quadrant N × N, starting at (0, 0) and ending at some point (n, k), using the same steps as before.Let J n,k (β, γ) be the generating polynomial for partial Motzkin paths from (0, 0) to (n, k), in which each rise gets weight 1, each level step at height i gets weight γ i , and each fall from height i gets weight β i .We therefore have an infinite unit-lower-triangular array J = J n,k (β, γ) n,k≥0 in which the first (k = 0) column displays the ordinary Jacobi-Rogers polynomials J n,0 = J n .It is immediate from the definition of J n,k that the matrix J is the output matrix O(Π) corresponding to the tridiagonal production matrix which arises from splitting a Motzkin path of length n+n into its first n steps and its last n steps, and then imagining the second part run backwards: the factor k i=1 β i arises from the fact that when we reversed the path we interchanged rises with falls and thus lost a factor k i=1 β i for those falls that were not paired with rises.The identity (2.29) can be written in matrix form as (2.27).
Remarks. 1.The reversal argument employed in this proof can be rewritten purely algebraically as follows: Note first that the tridiagonal matrix (2.25) satisfies where D = diag(1, β 1 , β 1 β 2 , . ..). (Here we work in the ring Z[β, β −1 , γ] of Laurent polynomials in β.) On the other hand, it is a general fact [65] that if Π is a production matrix and O 0 (Π) is the zeroth-column sequence of O(Π), then And finally, it is a general fact [65] that if M is an invertible lower-triangular matrix satisfying Putting all this together, we have as asserted in Proposition 2.2.Obviously, this proof relies crucially on the fact that the reversal of a Motzkin path is again a Motzkin path, or equivalently on the fact that the production matrix Π is tridiagonal and therefore symmetric up to a diagonal similarity transformation (2.30).
Taking the determinant of the n × n leading principal submatrix on both sides of (2.27), we obtain a classical formula [82,Theorem 51.1] for the Hankel determinants of a J-fraction: In particular, if R = R and all β i > 0, then the Hankel matrix H ∞ (J ) is positivedefinite, which implies that the underlying sequence (J n (β, γ)) n≥0 is a Hamburger moment sequence with a representing measure of infinite support [1,3,60,63].

Branched continued fractions
In this subsection we give a very brief introduction to the theory of branched continued fractions, limiting attention for simplicity to branched S-fractions; our treatment follows [56], where many more details and applications can be found.
Fix an integer m ≥ 1.An m-Dyck path [7,12,56,57] is a path in the upper half-plane Z × N, starting and ending on the horizontal axis, using steps (1, 1) ["rise" or "up step"] and (1, −m) ["m-fall" or "down step"].More generally, an m-Dyck path at level k is a path in Z × N ≥k , starting and ending at height k, using steps (1, 1) and (1, −m).Since the number of up steps must equal m times the number of down steps, the length of an m-Dyck path must be a multiple of m + 1.
Now let α = (α i ) i≥m be an infinite set of indeterminates.Then [56] the m-Stieltjes-Rogers polynomial of order n, denoted S Let n (α) t n be the ordinary generating function for m-Dyck paths with these weights; and more generally, let f k (t) be the ordinary generating function for m-Dyck paths at level k with these same weights.(Obviously f k is just f 0 with each α i replaced by α i+k ; but we shall not explicitly use this fact.)Then straightforward combinatorial arguments [56, Section 2.3], similar to those used in the proof of Theorem 2.1, lead to the functional equation or equivalently Iterating (2.35), we see immediately that f k is given by the branched continued fraction and in particular that f 0 is given by the specialization of (2.36) to k = 0. We shall call the right-hand side of (2.36) an m-branched Stieltjes-type continued fraction, or m-branched S-fraction for short.
Remark.In truth, we hardly ever use the branched continued fraction (2.36); instead, we work directly with the m-Dyck paths and/or with the recurrence (2.34)/(2.35) that their generating functions satisfy.
We now generalize these definitions as follows.A partial m-Dyck path is a path in the upper half-plane Z × N, starting on the horizontal axis but ending anywhere, using steps (1, 1) ["rise"] and ( • L(s 1 , s 2 , . ..) is the lower-bidiagonal matrix with 1 on the diagonal and s 1 , s 2 , . . . on the subdiagonal: • U (s 1 , s 2 , . ..) is the upper-bidiagonal matrix with 1 on the superdiagonal and s 1 , s 2 , . . . on the diagonal: Then the production matrix for the triangle S (m) is that is, the product of m factors L and one factor U [56, Proposition 8.2].
Let us remark, finally, that there is (as far as I know) no analogue of Proposition 2.2 for m-S-fractions with m > 1, since the reversal of an m-Dyck path is not an m-Dyck path.
3 Production matrix for a sequence of monic polynomials In this section we prove some very elementary -but important -relations between sequences of monic polynomials, the linear recurrences they satisfy, and their dual sequences of linear functionals.All of these properties will be re-expressed in a convenient matrix form, using the theory of production matrices (Section 2.2).We conclude this section with some more delicate matters concerning "expectation values" of products of polynomials, culminating in Open Problem 3.10.

Linear functionals
Let R[x] be the ring of polynomials in one indeterminate x, with coefficients in To each linear functional L : R[x] → R there is naturally associated a sequence ( n ) n≥0 of elements of R, defined by n = L(x n ).And conversely, to every sequence ( n ) n≥0 there is associated a unique linear functional We call ( n ) n≥0 the moment sequence of the linear functional L. Now let (L k ) k≥0 be a sequence of such linear functionals.We form the matrix A = (a nk ) n,k≥0 whose columns are the moment sequences of these linear functionals, i.e. a nk = L k (x n ).We call A the moment matrix for the sequence (L k ) k≥0 of linear functionals.We say that the sequence (L k ) k≥0 is normalized in case the matrix A is unit-lower-triangular, i.e.L k (x n ) = 0 for n < k and L k (x k ) = 1.

Sequences of monic polynomials
By a sequence of monic polynomials we mean a sequence (P n (x)) n≥0 of polynomials (with coefficients in R) such that P n (x) has degree n and leading coefficient 1.We can assemble the coefficients of these polynomials into a unit-lower-triangular matrix B = (b nk ) n,k≥0 by writing There is obviously a one-toone correspondence between unit-lower-triangular matrices and sequences of monic polynomials.We call B the coefficient matrix for the sequence (P n (x)) n≥0 of polynomials.
Now let A = (a nk ) n,k≥0 be the inverse matrix to B, i.e.A = B −1 .Then we obviously have a nk P k (x).

Duality
Let P = (P n (x)) n≥0 be a sequence of monic polynomials, and let L = (L k ) k≥0 be a sequence of linear functionals.We say that P and L are dual to each other in case L k (P n (x)) = δ kn for all k, n ≥ 0. The fundamental result concerning such duality is very simple: Proposition 3.1 (Sequence of monic polynomials and its dual sequence of linear functionals).Given any sequence (P n (x)) n≥0 of monic polynomials, there exists a unique sequence (L k ) k≥0 of linear functionals that satisfies L k (P n (x)) = δ kn , and it is normalized.
Conversely, given any normalized sequence (L k ) k≥0 of linear functionals, there exists a unique sequence (P n (x)) n≥0 of monic polynomials that satisfies L k (P n (x)) = δ kn .
The relation between these sequences is: The moment matrix A of the sequence (L k ) k≥0 and the coefficient matrix B of the sequence (P n (x)) n≥0 are inverses of each other.
Proof.Using P n (x) = n j=0 b nj x j and L k (x j ) = a jk , we see that the condition L k (P n (x)) = δ kn is equivalent to the matrix equation BA = I.If (P n (x)) n≥0 is a sequence of monic polynomials, then B is unit-lower-triangular, and BA = I has the unique solution A = B −1 ; moreover, A is unit-lower-triangular.And conversely, if (L k ) k≥0 is a normalized sequence of linear functionals, then A is unit-lower-triangular, and BA = I has the unique solution B = A −1 ; moreover, B is unit-lower-triangular.

Linear recurrence ←→ production matrix
The next result, which is only slightly more complicated, connects a sequence (P n (x)) n≥0 of monic polynomials with the unique linear recurrence (of a certain standard form) that defines it: Proposition 3.2 (Sequence of monic polynomials and its defining recurrence).Given any sequence (P n (x)) n≥0 of monic polynomials, there exists a unique unit-lower-Hessenberg matrix Π = (π nk ) n,k≥0 such that or equivalently And conversely, given any unit-lower-Hessenberg matrix Π = (π nk ) n,k≥0 , there exists a unique sequence (P n (x)) n≥0 of polynomials satisfying (3.1)/(3.2) with the initial condition P 0 (x) = 1, and it is monic.
The relation between these objects is: The coefficient matrix B of the sequence Proof.Let (P n (x)) n≥0 be a sequence of monic polynomials with coefficient matrix B. Substituting P n (x) = n j=0 b nj x j into (3.2) and extracting the coefficient of x j , we see that ( 2. The recurrence (3.2) can also be written in vector form as xP = ΠP, where P is the column vector whose entries are the sequence (P n (x)) n≥0 .Iterating this, we see that x r P = Π r P for any integer r ≥ 0, or concretely It then follows that for any polynomial q(x) we have q(x) P n (x) = n+deg q j=0 (q(Π)) nj P j (x) . (3.6) The monic polynomials P n (x) defined by the recurrence (3.1) are also the characteristic polynomials of the leading principal submatrices of the production matrix Π.To state this result, let us introduce notation as follows: For any matrix A = (a ij ) i,j≥0 , we write A n = (a ij ) 0≤i,j≤n−1 for its n × n leading principal submatrix, and ∆ n (A) = det A n for its n × n leading principal minor, with the convention ∆ 0 (A) = 1.When the matrix is lower-Hessenberg, these leading principal minors satisfy a recurrence that is reasonably well known [32, pp.251-252] [84, pp.410-411], though perhaps not as well known as it should be:

Lemma 3.3 (Leading principal minors of a lower-Hessenberg matrix).
The leading principal minors ∆ n (H) of a lower-Hessenberg matrix H = (h ij ) i,j≥0 satisfy the recurrence (blank entries are zero), so that its determinant is For tridiagonal matrices, (3.7) becomes a three-term recurrence that is much better known [35, p. 35].
Proposition 3.4 (Monic polynomials as characteristic polynomials of production matrix).Let Π = (π nk ) n,k≥0 be a unit-lower-Hessenberg matrix, and let (P n (x)) n≥0 be the sequence of monic polynomials defined by the recurrence (3.1) with the initial condition P 0 (x) = 1.
Remarks. 3. When R = R or C, it follows from Proposition 3.4 that the zeros of P n (x) are the eigenvalues of Π n .This suggests that the asymptotic zero distribution of the polynomials P n (x) as n → ∞ should be related to the spectral properties of the infinite matrix Π acting on a suitable space of sequences (for instance, 2 (N)).See e.g.[6,39,40,58,87].
Finally, the recurrence (3.1) also leads to a combinatorial formula, due to Viennot [80, p. III-16], for the matrix elements of B = O(Π) −1 in terms of those of Π: Corollary 3.5 (Viennot [80]).Let Π = (π nk ) n,k≥0 be a unit-lower-Hessenberg matrix, and let B = (b nj ) n,j≥0 be given by B = O(Π) −1 .Then b nj is the sum over partitions of {0, 1, . . ., n − 1} into zero or more intervals [k, ] (k ≤ ) and exactly j empty sites, with a weight −π k for each interval [k, ]. 2Proof.Write b nj for the quantity defined in the Corollary, and define P n (x) = n j=0 b nj x j .Then P n (x) is the sum over partitions of {0, 1, . . ., n − 1} into zero or more intervals [k, ] (k ≤ ) and zero or more empty sites, with a weight −π k for each interval [k, ] and a weight x for each empty site.And it is easy to see, by considering the status of the vertex n in P n+1 (x), that the sequence ( P n (x)) n≥0 satisfies the same recurrence (3.1) as is satisfied by (P n (x)) n≥0 , with the same initial condition P 0 (x) = P 0 (x) = 1.So P n (x) = P n (x).

Summary
To summarize the results obtained thus far: There is a one-to-one correspondence between sequences (P n (x)) n≥0 of monic polynomials (with coefficient matrix B), their dual sequences (L k ) k≥0 of linear functionals (with moment matrix A), and their defining linear recurrences (3.1)/(3.2) (with production matrix Π); and these correspondences are given by (2.11).

Expectation values of products
Fix now a unit-lower-Hessenberg matrix Π: this defines (by Proposition 3.2) a sequence (P n (x)) n≥0 of monic polynomials with coefficient matrix B = O(Π) −1 , which in turn defines (by Proposition 3.1) a dual sequence (L k ) k≥0 of linear functionals with moment matrix A = O(Π).Our goal is to find a general combinatorial or algebraic formula for quantities of the form L k (q(x) P m (x) P n (x)), where q(x) is a polynomial, in terms of the coefficients Π.By analogy to probability theory, we refer colloquially to quantities like L k (q(x) P m (x) P n (x)) as "expectation values".
In the tridiagonal case with k = 0, Viennot [80, p. I-15, Proposition I.17] found a beautiful formula for these expectation values: Proposition 3.6 (Viennot [80]).When the unit-lower-Hessenberg matrix Π is tridiagonal, we have for any polynomial q(x).In particular, when q(x) = 1 we have the orthogonality relation L 0 (P m (x) P n (x)) = h n δ mn (3.11) with the normalizing constant Please note [80, p. I-15] that the right-hand side of (3.10) with q(x) = x r , namely π 10 π 21 • • • π n,n−1 (Π r ) mn , can be interpreted as the total weight for Motzkin paths of length r + m + n from height 0 → 0 in which the first m steps are "up" steps (getting weight 1) and the last n steps are "down" steps (getting weight π n,n−1 • • • π 10 ).Now, in a Motzkin path that starts and ends at the same height, each "up" step i → i + 1 can be paired with a "down" step i + 1 → i; it follows that the weight of such a path equals the weight of the reversed path.These considerations show that the right-hand side of (3.10) is indeed symmetric in m ↔ n.
Viennot [80, pp.I-16-I-19] proved Proposition 3.6 by a rather intricate combinatorial argument; here we give a simple algebraic proof: Proof of Proposition 3.6.Let = ( n ) n≥0 be the moment sequence of the linear functional L 0 , i.e. (3.12) And let H ∞ ( ) = ( i+j ) i,j≥0 be the Hankel matrix associated to the sequence .Since the unit-lower-Hessenberg matrix Π is tridiagonal, it is of the form (2.25) with γ n = π nn and β n = π n,n−1 .Proposition 2.2 therefore gives where or concretely But the left-hand side of (3.15) is exactly L 0 (P m (x) P n (x)).This proves (3.11).
Hélder Lima [45] has pointed out to me that Proposition 3.6 can be extended to k = 0 as follows: Proposition 3.7.When the unit-lower-Hessenberg matrix Π is tridiagonal, we have and hence L k (q(x) P m (x) for any polynomial q(x), where Proof.Put m = 0 in (3.10), rename n as k, and take q(x) = x n : this gives for any polynomial q(x).This proves (3.18).Then replace q(x) by q(x) P m (x) P n (x) and use (3.10); this proves (3.19).
Remark.The meaning of h −1 k in the formulae (3.18) and (3.19) requires some clarification.The formulae obviously hold if h k is invertible in the ring R. In particular this is the case if we consider Π = {π ij } i≥j≥0 to be indeterminates and we work in the ring Z[Π, Π −1 ] of Laurent polynomials.Now suppose that k ≤ n: then both sides of (3.19) are in fact polynomials in Π, so the identity holds when the π ij are specialized to arbitrary elements in an arbitrary commutative ring.I am not sure what happens when k > n.
We would now like to generalize these results to the non-tridiagonal case.(Viennot [80, alludes to this as an open problem.)For the moment I have only the following result: Theorem 3.8.With (P n (x)) n≥0 and (L k ) k≥0 defined as above, we have and more generally where H k is the -shifted Hankel matrix of the k th column of A: Proof.Since (P n (x)) n≥0 and (L k ) k≥0 are a dual pair, we have L k (P n (x)) = δ kn .
Applying L k to (3.5), we obtain L k (x P n (x)) = (Π ) nk ; this proves (3.22).Now use P m (x) = m j=0 b mj x j ; inserting this into (3.22)gives (3.23a).On the other hand, we have L k (x n ) = a nk .Inserting the representing equation in terms of B for both P m (x) and P n (x), we have k is defined by (3.24).This proves (3.23b).
Remark.The identity (3.22) is contained in the thesis of Roblet [59, p. 153, Proposition 78] in the special case where Π is (d, 1)-banded and 0 ≤ k ≤ d − 1 (but this is no real loss of generality, since we can take d → ∞).Roblet's proof was combinatorial, following the model of Viennot [80, pp. I-15-I.19].
It is an immediate consequence of Theorem 3.8 that certain matrix elements have to vanish: Corollary 3.9.In the situation of Theorem 3.8: Proof.(a) is an immediate consequence of (3.23a) together with the facts that B is lower-triangular and Π is lower-Hessenberg.(b) The vanishing for k < n − d( + m) is likewise an immediate consequence of (3.23a) together with the facts that B is lower-triangular and Π is (d, 1)-banded.The vanishing for k < m − d( + n) then follows from the symmetry m ↔ n of the left-hand side of (3.23a).
In particular, if Π is (d, 1)-banded, then L k (x P n (x)) = 0 whenever n > d +k.So the sequence (P n (x)) n≥0 is d-orthogonal with respect to the sequence L 0 , . . ., L d−1 of linear forms, in the sense defined in the Introduction.But in fact this vanishing for n > d + k holds for all k ≥ 0, not just for k ≤ d − 1.
Unfortunately equation (3.23a) is not very nice, because it is not manifestly symmetric in m ↔ n; moreover, it is not expressed solely in terms of Π.The equation (3.23b) has the m ↔ n symmetry, but it is still not expressed solely in terms of Π. Combining (3.23a) with Corollary 3.5 gives an explicit formula for L k (x P m (x) P n (x)) in terms of the matrix elements of Π; but this formula is rather complicated, and it also fails to make manifest the symmetry m ↔ n.Combining (3.23b) and (3.24) with Corollary 3.5 gives an even more complicated formula for L k (x P m (x) P n (x)) in terms of the matrix elements of Π, which at least is manifestly symmetric in m ↔ n.But none of these formulae seem really satisfactory.What we really want is something that looks more like Propositions 3.6 and 3.7 and that reduces to them when Π is tridiagonal.My hope is that there might be some cancellations between terms in the expansion of (3.23a), so that the final result can perhaps be stated in a simpler way.I therefore conclude by stating the main unsolved problem of this paper: Open Problem 3.10.Find a more satisfactory formula for L k (x P m (x) P n (x))ideally one that resembles Propositions 3.6 and 3.7 and that reduces to them when Π is tridiagonal.Even the special case k = = 0 would be of great interest.

Sequence of monic polynomials orthogonal to a sequence of linear functionals
Let Γ = (Γ k ) k≥0 be a sequence of linear functionals on R[x], with moment matrix Γ = (γ nk ) n,k≥0 given by γ nk = Γ k (x n ).And let P = (P n (x)) n≥0 be a sequence of monic polynomials, with unit-lower-triangular coefficient matrix B = (b nj ) n,j≥0 given by P n (x) = n j=0 b nj x j .We say that P is orthogonal to Γ in case each we see that P is orthogonal to Γ if and only if BΓ vanishes below the diagonal, or in other words BΓ is an upper-triangular matrix U , or equivalently Γ = B −1 U .We record this simple fact: Proposition 4.1 (Orthogonality between a sequence of linear functionals and a sequence of monic polynomials).Let Γ = (Γ k ) k≥0 be a sequence of linear functionals on R[x], with moment matrix Γ; and let P = (P n (x)) n≥0 be a sequence of monic polynomials, with unit-lower-triangular coefficient matrix B. Then P is orthogonal to Γ if and only if the matrix BΓ is upper-triangular, or equivalently if there exists an upper-triangular matrix U such that Γ = B −1 U .
We can interpret this result in two ways, depending on whether we start from Γ or from P: Starting from Γ. Suppose that the moment matrix Γ has a factorization Γ = LU where L is unit-lower-triangular and U is upper-triangular.Then B = L −1 is the coefficient matrix for a sequence P of monic polynomials that is orthogonal to Γ. Furthermore, the sequence L = (L k ) k≥0 of linear functionals dual to P (given by Proposition 3.1) then has moment matrix A = B −1 = L.It follows that Γ k = k j=0 (U T ) kj L j ; and if the diagonal elements of U are invertible, then L k = k j=0 (U −T ) kj Γ j .In other words, each Γ k is a linear combination of L 0 , . . ., L k ; and if the diagonal elements of U are invertible, then each L k is a special linear combination of Γ 0 , . . ., Γ k , namely, one that makes the first k elements of its moment sequence zero and the next element 1.
Starting from P. Let P = (P n (x)) n≥0 be a sequence of monic polynomials, with coefficient matrix B. Then there is a canonically associated sequence of linear functionals with respect to which P is orthogonal, namely, the dual sequence L = (L k ) k≥0 , with moment matrix A = B −1 .But P is also orthogonal with respect to any sequence Γ = (Γ k ) k≥0 of linear functionals whose moment matrix Γ is of the form Γ = AU , where U is any upper-triangular matrix; or in other words, each Γ k is an arbitrary linear combination of L 0 , . . ., L k .
Remark.When R = R, it is often the case in applications that L 0 is the moment functional of a positive measure.But L k for k ≥ 1 cannot be the moment functional of a positive measure, because the zeroth component of its moment sequence vanishes but the sequence is not identically zero (since the kth component is 1).Analyst readers may be interested in the following problem: Given a unit-lower-triangular matrix A whose zeroth column is the moment sequence of a positive measure, find an upper-triangular matrix U such that every column of Γ = AU is the moment sequence of a positive measure; or in other words, find linear combinations Γ k of L 0 , . . ., L k that are all moment functionals of positive measures.This problem of course has trivial solutions: we could take all the Γ k to be zero, or to be equal to L 0 .But suppose we further insist that all the functionals Γ 0 , Γ 1 , . . .be linearly independent.Then the problem seems to be nontrivial.
Since Γ k is a linear combination of L 0 , . . ., L k whenever P is orthogonal to Γ, it also follows that the result of Corollary 3.9(b) holds with L k replaced by Γ k : Corollary 4.2.Let Π = (π nk ) n,k≥0 be a unit-lower-Hessenberg matrix, let P = (P n (x)) n≥0 be the sequence of monic polynomials defined by (3.1)/(3.2) with initial condition P 0 (x) = 1, and let Γ = (Γ k ) k≥0 be any sequence of linear functionals such that P is orthogonal to Γ.
In view of the key role played here by the factorization Γ = LU when we start from Γ, it is now appropriate to recall some simple facts concerning the existence and uniqueness of LU factorizations for matrices over a commutative ring [66].Let us say that a square matrix Γ has a weak LU factorization if Γ = LU where L is lower-triangular and U is upper-triangular, and an LU factorization if Γ = LU where L is unit-lower-triangular and U is upper-triangular.We then have [66]: Proposition 4.3 (LU factorization for matrices over a commutative ring).Let Γ be an n × n matrix with entries in a commutative ring R, and let ∆ 1 , . . ., ∆ n be its leading principal minors.
(a) If Γ has a weak LU factorization, then we must have Since [66] is not yet publicly available, we include a proof of Proposition 4.3 in the Appendix.
Taking n → ∞ in Proposition 4.3(d,e) and applying it to the situation considered in Proposition 4.1, we conclude: Corollary 4.4 (Existence and uniqueness of a sequence of monic polynomials orthogonal to a given sequence of linear functionals).Let R be a commutative ring, and let Γ = (Γ k ) k≥0 be a sequence of linear functionals on R[x], with moment matrix Γ.
(a) If none of the leading principal minors ∆ 1 , ∆ 2 , . . . of Γ is zero or a divisor of zero, then there is at most one sequence of monic polynomials orthogonal to Γ. (In particular this holds if R is an integral domain and ∆ 1 , ∆ 2 , . . .= 0.) (b) If all of the leading principal minors ∆ 1 , ∆ 2 , . . . of Γ are invertible in R, then there is exactly one sequence of monic polynomials orthogonal to Γ. (In particular this holds if R is a field and ∆ 1 , ∆ 2 , . . .= 0.)

Application to ordinary orthogonal polynomials
Let us begin by showing how the general theory from the preceding section applies to ordinary orthogonal polynomials.
Fix a linear functional L, with moment sequence = ( n ) n≥0 given by n = L(x n ).And let us choose Γ k to be the k-shift of L: that is, Γ k (x n ) def = L(x n+k ).Then the moment matrix Γ of the sequence Γ = (Γ k ) k≥0 of linear functionals is the Hankel matrix H ∞ ( ) = ( i+j ) i,j≥0 associated to the sequence : that is, γ nk = n+k .And a sequence P = (P n (x)) n≥0 of monic polynomials is orthogonal to Γ in case L(x k P n (x)) = 0 for 0 ≤ k ≤ n − 1, i.e. precisely when P is a sequence of monic orthogonal polynomials in the usual sense associated to the linear functional L [14, Chapter 1].In particular, by Corollary 4.4(b), such a sequence P exists (and is unique) whenever R is a field and all the leading principal minors ∆ 1 , ∆ 2 , . . . of Γ are nonzero.
Let us now relate this to production matrices and classical continued fractions.It is known [82,Theorem 51.1] [80, p. IV-17, Corollaire 7 and p. V-5, Proposition 1] that if R is a field and all the leading principal minors ∆ 1 , ∆ 2 , . . . of the Hankel matrix Γ = H ∞ ( ) are nonzero, then there exists a classical J-fraction that represents the ordinary generating function of the sequence that underlies this Hankel matrix, i.e.
in the sense of formal power series, with coefficients γ 0 , γ 1 , . . .∈ R and β 1 , β 2 , . . .∈ R \ {0}.In fact, the J-fraction coefficients are connected to the leading principal minors by where ∆ −1 def = 1 [this follows from (2.33)], while the γ n are given by other determinants involving the moments [80, Sections IV.3 and V.1]. 3his J-fraction has all the properties described in Section 2. But it is pleasing to see all these classical facts come together as consequences of a simple algebraic theory.

Application to multiple orthogonal polynomials
Let us now apply the general theory from the Section 4 to the multiple orthogonal polynomials of type II along an increasing nearest-neighbor path in N r .
If L is a linear functional on R[x] and k is a nonnegative integer, we denote by L k the k-shift of L: that is, L k (x n ) def = L(x n+k ).Now fix an integer r ≥ 1, and fix linear functionals L (1) , . . ., L (r) on R[x].(In the analytical setting, we will have R = R, and L (1) , . . ., L (r) will be the moment functionals associated to the positive measures µ 1 , . . ., µ r .)And let (P n (x)) n∈N r be the multiple orthogonal polynomials of type II associated to the linear functionals L (1) , . . ., L (r) , which we here assume to exist.Now let j 1 , j 2 , . . .be an infinite sequence of elements of {1, . . ., r}, and define a sequence (n k ) k≥0 of multi-indices in N r by n k = k i=1 e j i .They satisfy |n k | = k and describe an increasing nearest-neighbor path in N r in which the ith step is along direction j i .Now let P k (x) def = P n k (x) be the multiple orthogonal polynomial of type II along this path in N r .Let m i = (n i−1 ) j i = (n i ) j i − 1 be the number of indices j 1 , . . ., j i−1 that equal j i .Then P k (x) is orthogonal to the linear functionals L 1 , . . ., L k , where the "new" linear functional appearing at stage k is i.e.
Now set Γ k = L ,k+1 : we then see that the sequence P = ( P k (x)) k≥0 is orthogonal to the sequence Γ = (Γ k ) k≥0 in the sense of the preceding section.
In the case r = 2, the formulae for L 0 and L 1 (but not the higher L k ) can be found already in the thesis of Drake [27, Theorem 2.4.1].
7 Some examples 7.1 Bessel K ν weights ⇒ Rising-factorial moments Two decades ago, Van Assche and Yakubovich [76] studied the multiple orthogonal polynomials of types I and II associated to a pair of measures (that is, r = 2) in which the weights are modified Bessel functions of the second kind, K ν (x) [83, p. 78], multiplied by powers of x.We shall follow their paper closely, but change the notation to make the formulae more symmetrical.
In fact, we can go farther and compute the full output matrix O(Π), i.e. compute the generalized 2-Stieltjes-Rogers polynomials S n,k (α) for the coefficients α given by (7.10).This was not done in [56], but we can do it here: Proposition 7.1 (Generalized 2-Stieltjes-Rogers polynomials associated to the rising-factorial moments).The output matrix O(Π) = S (2) (α) corresponding to the production matrix (7.5) is Proof.Let a nk be the right-hand side of (7.12); we need to show that it satisfies the recurrence (2.8) when the π ik are given by (7.5).That is, we need to show that This is a tedious but straightforward computation: it is convenient to pull out from the right-hand side a factor (n−1)! (n−k)!(k+2)!(a 1 + k + 2) n−k−3 (a 2 + k + 2) n−k−3 and then evaluate the remaining polynomial expression.
4. More generally, one can consider cases in which the moments of the measures µ 1 , . . ., µ r are ratios of products of rising factorials, with p factors in the numerator and q factors in the denominator.Then the ordinary generating function of the moments of µ 1 is a hypergeometric function F p+1 q with a p+1 = 1, and the recurrence relation for the stepline polynomials can be compared with the branched continued fractions in [56,Theorems 14.3,14.5 and 14.6].In all these branched continued fractions, r = max(p, q).Lima and Loureiro have considered the cases (p, q) = (2, 1) [47] and (p, q) = (2, 2) [48]; and Lima [46] has very recently considered the case of general (p, q).This quadridiagonal production matrix Π plays a central role in our forthcoming work [21] on the coefficientwise Hankel-total positivity of the Laguerre polynomials.When α = −1 (Lah polynomials) it arises from a 2-branched S-fraction, as found already in [55]; when α = 0 (rook polynomials) it arises from a modified 2-branched S-fraction (see [21]).
Coussement and Van Assche [18] also gave an explicit formula for the multiple orthogonal polynomials of type II along the stepline [18, Theorem 10 and Corollary 2].After translating from their notation to ours, it is It is curious that Laguerre polynomials occur here too.

Final remarks
I suspect that the foregoing examples are just the tip of the iceberg, and that the connection between multiple orthogonal polynomials, production matrices and branched continued fractions will be fruitful in both directions.For instance, using known techniques (such as vector Pearson equations [9,17,18,26,47,48,76]) it may be possible to devise new examples of multiple orthogonal polynomials; these will then automatically provide a production matrix for the sequence of moments of µ 1 ; this production matrix will in turn automatically arise from a branched J-fraction [56,, and in some cases this branched J-fraction may arise from contraction of a branched S-fraction [56,Propositions 7.2 and 7.6].And conversely, combinatorial or algebraic methods leading to new production matrices or branched continued fractions may point the way to new examples of multiple orthogonal polynomials.
(m) n (α), is the generating polynomial for m-Dyck paths of length (m + 1)n in which each rise gets weight 1 and each m-fall from height i gets weight α i .Clearly S (m) n (α) is a homogeneous polynomial of degree n with nonnegative integer coefficients.
1, −m) ["m-fall"].A partial m-Dyck path starting at (0, 0) must stay always within the set V m = {(x, y) ∈ Z × N : x = y mod m + 1}.Now let α = (α i ) i≥m be an infinite set of indeterminates, and let S (m) n,k (α) be the generating polynomial for partial m-Dyck paths from (0, 0) to ((m + 1)n, (m + 1)k) in which each rise gets weight 1 and each m-fall from height i gets weight α i .We call the S (m) n,k the generalized m-Stieltjes-Rogers polynomials.Obviously S (m) n,k is nonvanishing only for 0 ≤ k ≤ n, and S (m) n,n = 1.We therefore have an infinite unit-lower-triangular array S (m) = S (m) n,k (α) n,k≥0 in which the first (k = 0) column displays the ordinary m-Stieltjes-Rogers polynomials S (m) n,0 = S (m) n .The production matrix for the triangle S (m) was found in [56, sections 7.1 and 8.2].We begin by defining some special matrices M = (m ij ) i,j≥0 : where a | b denotes that a divides b in the ring R. (b) If Γ has a weak LU factorization in which none of the diagonal elements of L or U is a zero or a divisor of zero, then none of ∆ 1 , . . ., ∆ n is zero or a divisor of zero.(c) If Γ has a weak LU factorization in which all of the diagonal elements of L and U are invertible in R, then ∆ 1 , . . ., ∆ n are invertible in R. Conversely, (d) If none of ∆ 1 , . . ., ∆ n−1 is a zero or a divisor of zero, then Γ has at most one LU factorization.(In particular this holds if R is an integral domain and ∆ 1 , . . ., ∆ n−1 = 0.) (e) If ∆ 1 , . . ., ∆ n−1 are invertible in R, then Γ has exactly one LU factorization.(In particular this holds if R is a field and ∆ 1 , . . ., ∆ n−1 = 0.)