Properties of matrix orthogonal polynomials via their Riemann-Hilbert characterization

Symmetry, Integrability and Geometry: Methods and Applications SIGMA 7(2011), 098, 31 pages Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization F. Alberto GR ¨ UNBAUM †, Manuel D. DE LA IGLESIA ‡ and Andrei MART´ INEZ-FINKELSHTEIN § †Department of Mathematics, University of California, Berkeley, Berkeley, CA 94720 USA E-mail: grunb[email protected]erkeley.edu URL: http://math.berkeley.edu/~grunbaum/ ‡Departamento de An´alisis Matem´atico, Universidad de Sevilla, Apdo (P. O. BOX) 1160, 41080 Sevilla, Spain E-mail: [email protected] URL: http://euler.us.es/~mdi29/ §Departamento de Estad´ıstica y Matem´atica Aplicada, Universidad de Almer´ıa, 04120 Almer´ıa, Spain E-mail: andr[email protected] URL: http://www.ual.es/~andrei/ Received June 09, 2011, in final form October 20, 2011; Published online October 25, 2011 http://dx.doi.org/10.3842/SIGMA.2011.098 Abstract. We give a Riemann–Hilbert approach to the theory of matrix orthogonal polynomials. We will focus on the algebraic aspects of the problem, obtaining difference and differential relations satisfied by the corresponding orthogonal polynomials. We will show that in the matrix case there is some extra freedom that allows us to obtain a family of ladder operators, some of them of 0-th order, something that is not possible in the scalar case. The combination of the ladder operators will lead to a family of second-order differential equations satisfied by the orthogonal polynomials, some of them of 0-th and first order, something also impossible in the scalar setting. This shows that the differential properties in the matrix case are much more complicated than in the scalar situation. We will study several examples given in the last years as well as others not considered so far. Key words: matrix orthogonal polynomials; Riemann–Hilbert problems 2010 Mathematics Subject Classification: 42C05; 35Q15 1 Introduction The theory of matrix orthogonal polynomials on the real line (MOPRL) has its foundations in the seminal papers of Krein [42,43] (see also their account in the book of Berezans’ki˘ı [2]). For further historical background and analytic results, the reader is referred to the survey [11] and to Chapter 4 of [51]. In many aspects the MOPRL resemble their scalar counterparts, especially where the proofs are based on basic properties of Hilbert spaces. Nevertheless, the non-commutativity of matrix multiplication and the existence of non-zero singular matrices add features to the theory that make MOPRL an interesting object of study. Moreover, many problems for scalar polynomials are better understood or recast in terms of some matrix polynomials, see, for instance, [26]. Recall that a matrix polynomial of degree ≤nin CN×Nand a scalar variable xcan be defined as an expression of the form Anxn+· · · +A1x+A0, 2 F.A. Gr¨unbaum, M.D. de la Iglesia and A. Mart´ınez-Finkelshtein where Aj’s are constant matrices in CN×N. In what follows we consider Nfixed, and denote by Pnthe family of matrix polynomials of degree ≤nin CN×N, as well as by P:= Sn≥0Pn. We use preferably boldface letters to denote matrices, and standard font for scalars. We also use INfor the N×Nidentity matrix, omitting the explicit reference to its dimension when it cannot lead anyone into confusion. In this paper we consider some aspects of the MOPRL theory, assuming that the orthogonality is given by an absolutely continuous measure on the line. More precisely, our starting point is a weight W= (Wij) : (a, b)→GL(N, R), defined and positive definite on a finite or infinite interval (a, b)⊂R. We will assume that all Wij and W0 ij have finite moments: Zb a |x|nW(x)dx < ∞,Zb a |x|nW0(x)dx < ∞, n ∈N0:= N∪ {0}, where the integration of a matrix function is applied entry-wise. For any two P,Q∈P, the weight Winduces two matrix-valued “inner products”, (P,Q)W=Zb a P(x)W(x)Q∗(x)dx, and hP,QiW=Zb a Q∗(x)W(x)P(x)dx = (Q∗,P∗)W, where the asterisk denotes the conjugate transpose (or Hermitian conjugate) of a matrix. Due to this connection between both inner products, we restrict our attention to (·,·)W. We define also the norm kPkW:= (Tr hP,PiW)1/2, and assume that Wis non-trivial, in the sense that kPkW>0 for every non-zero matrix polynomial P. In this case (see [11, Lemma 2.3] or [51, Proposition 4.2.3]), (P,P)Wis nonsingular for every non-zero polynomial P, and we can easily implement a matrix analogue of the Gram–Schmidt orthogonalization procedure, which yields a unique sequence ( b Pn)nof monic orthogonal polynomials such that b P0=I, b Pn(x) = xnI+ n−1 X j=0 an,jxj,(b Pn,Q)W=0for every Q∈Pn−1, n ∈N,(1.1) as well as matrix polynomials (Pn)n, “orthonormal” with respect to W, such that Pn(x) = κnb Pn,(Pn,Pm)W=δn,mI.(1.2) Obviously, Pn’s are determined up to a unitary left factor, so we can speak about an equivalence class of orthonormal MOPRL, corresponding to the weight W. Additionally, if Bis a constant non-singular matrix, then f W(x) = BW (x)B∗ is also a weight, and ( e Pn)n= (PnB−1)nis the corresponding sequence of orthonormal polynomials. Hence, more than a single weight we consider an equivalence class of weights given by its representative W. In particular, without loss of generality we can assume that W(x0) = I at a point x0∈(a, b). If this class contains a diagonal matrix-valued function, we say that W Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization 3 reduces to scalar weights. A characterization of this fact, as shown in [23], is the commutativity condition W(x)W(y) = W(y)W(x), x, y ∈[a, b]. This situation is considered trivial, and usually omitted from consideration. Work in the last few years has revealed a number of explicit families of MOPRL; in many cases they are joint eigenfunctions of some fixed differential operator with matrix coefficients independent of the degree nof the MOPRL. This study was initiated in [21], but nontrivial examples had to wait until [23,29,34,36]. A solution to the classification problem for those W whose MOPRL are common eigenfunctions of some fixed differential operator remains elusive. There are by now two methods that have yielded nontrivial examples: a) the connection with matrix-valued spherical functions for symmetric spaces [34,36], and b) a combination of classical methods and some Lie algebra tools [23,24]. A necessary condition in terms of moments is given in [23], and a necessary and sufficient condition in terms of the “ad-conditions” is given in [38]. In this approach one gets an explicit formula for one differential operator. Concerning applications, the recurrent relations have been used in the study of certain quasi-birth-and-death processes [30,31,32,33], in which case the differential operator plays no role. Recently, new applications of these processes have been related to urn or Young diagrams models [37]. A big boost in the research of the scalar OPRL has been their Riemann–Hilbert (RH) characterization, introduced in the seminal paper of Fokas, Its and Kitaev [28], and complemented with a non-linear steepest descent analysis in a series of works of Deift, Zhou and collaborators [12,13,16,17]. This combination has allowed to establish extremely strong asymptotic results with applications in random matrix theory, in particular for unitary invariant ensembles [12,14], determinantal point processes [8,9,44,45], orthogonal Laurent polynomials [46,47], Painlev´e trascendents [10,41], the Toda lattices [15], to mention a few. The RH characterization, even without the steepest descent analysis, allows one to prove other algebraic and analytic properties of orthogonal polynomials, as it was illustrated in [13] for OP on the unit circle, and in [40, Chapter 22], for classical families of OPRL. One of the goals of this paper is to extend these considerations to MOPRL, something that, surprisingly enough, has not been explored in depth so far. We will show that the RH formulation, which has a very natural block-wise generalization to the MOPRL case, allows one to reveal some subtleties hidden in the matrix case. For instance, we will derive identities and difference-differential equations that were unknown even for some of the explicit families mentioned above. The RH technique can be also a natural method for obtaining the differential relations for MOPRL with respect to general weight matrices. There have been several papers considering the block-wise RH problem of the flavor similar to the one discussed in this paper, but for singular (namely, rank 1) weight matrices. For these weights the MOPRL are connected with a class of multiple orthogonal polynomials that find applications in the analysis of determinantal point processes and non-intersecting stochastic paths (see e.g. the works of Kuijlaars and collaborators [8,18,19,45]) or of the multicomponent 2D Toda lattice hierarchy [1]. This paper deals with strictly algebraic consequences of the RH formulation of matrix polynomials orthogonal with respect to a weight matrix supported on an interval [a, b]⊂R; the asymptotic analysis via the Deift–Zhou non-linear steepest descent method will not be discussed here. In Section 2we will discuss two dual RHPs, which are uniquely solved in terms of the MOPRL. Standard arguments related to RHP allow us to obtain the three-term recurrence relation and the differential identities (ladder operators) for MOPRL. Both yield also the so-called Lax pair for the polynomials, which is an over-determined system (2.28), whose compatibility conditions 4 F.A. Gr¨unbaum, M.D. de la Iglesia and A. Mart´ınez-Finkelshtein render further properties of MOPRL. All the results in this section have been previously obtained using different approaches, but the RH method will provide new and straightforward proofs. In Section 3we will consider a special transformation of the original RHP, based on factorizations of the weight matrix of the form W=T T ∗. The main goal is to obtain a RHP with constant jumps (independent of z), in order to simplify the differential relation obtained earlier. For simplicity, we will focus on the case where the support of the weight matrix is R. The extension of the method to weights supported on semi-infinite or finite intervals is simple, see a short discussion at the end of this paper. We exploit also the non-uniqueness of the weight factorization, which in the matrix case gives us some extra freedom and yields non-trivial relations. In particular, a family of ladder operators, some of them of 0-th order, is obtained, a phenomenon that is not possible in the scalar situation. Furthermore, by combining appropriately the family of ladder operators, we will get the 0-th, first and second order differential equations satisfied by the MOPRL. Some examples of lower order differential operators appear in [5,6,35]. The considerations so far have been completely general, applicable to any family of MOPRL. In Section 4we narrow the analysis of Sections 2and 3to some relevant examples of MOPRL supported on R, obtaining new results even in the cases studied previously in the literature. In the final stage of preparation of this manuscript we learned about the preprint [4] which also uses the Riemann–Hilbert approach for the analysis of the matrix orthogonal polynomial. Although there is some overlapping between the results contained in our Sections 2,3and in [4, Section 2], the focus of both contributions is different, and in this sense, complementary. 2 The Riemann–Hilbert problem for MOPRL 2.1 Formulation and basic properties In this section we discuss the Riemann–Hilbert problem (RHP) related to MOPRL with respect to a N×Nweight matrix Wsupported on an interval [a , b ]⊂R; this interval can be either bounded or unbounded. We assume Wcontinuous and non-vanishing on (a, b), and that at any finite endpoint of the support the weight Whas at worse a power-type singularity, that is, W is of the form W(z) = |z−c|γcf W(z), γc>−1,(2.1) where c∈ {a, b},c6=±∞, and f Wis a bounded, continuous and non-vanishing at z=c. This class of weights comprises all the examples considered so far in the literature. As a convention, in what follows we write the 2N×2Nmatrices partitioned into N×N blocks, as in (2.2) below. Recall that A∗stands for the Hermitian conjugate of the matrix A, as well as A−∗ = (A∗)−1. We adopt the convention that for any matrix-valued function P(z) of a complex variable z,P∗denotes the matrix-valued function obtained as P∗(z) := (P(¯z))∗. The RHP for MOPRL with respect to a weight matrix Wconsists in finding a matrix function Yn:C→C2N×2Nsuch that (Y1) Ynis analytic in C\[a, b]. (Y2) Ynhas on (a, b) continuous boundary values Yn +(resp., Yn −) from the upper (resp., lower) half plane, such that Yn +(x) = Yn −(x)INW(x) 0IN, x ∈(a, b).(2.2) Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization 5 (Y3) As z→ ∞, for every m∈Nwe have Yn(z) = I2N+ m X i=1 Yn i zi+O1/zm+1! znIN0 0z−nIN!(2.3) (where the asymptotic term O(1/zm+1) depends on n). (Y4) As z→c,c∈ {a, b},c6=±∞, we have Yn(z) = O(1) O(h(z)) O(1) O(h(z)) !, where h(z) =      |z−c|γc,if −1< γc<0, log |z−c|,if γc= 0, 1,if γc>0. (2.4) Remark 2.1. Usually the asymptotic condition at infinity (2.3) is stated for m= 1; however, it can be proved that (2.3) for m= 1 implies that this condition holds for every m∈N. Along with the RH problem (Y1)–(Y4) we can consider the following dual problem: finding a matrix function yn:C→C2N×2Nsuch that (y1) ynis analytic in C\[a, b]. (y2) yn +(x) = IN−W(x) 0INyn −(x) when x∈(a, b). (y3) As z→ ∞, for every m∈Nwe have yn(z) = z−nIN0 0znIN! I2N+ m X i=1 e Yn i zi+O1/zm+1!.(2.5) (y4) As z→c,c∈ {a, b},c6=±∞, we have yn(z) = O(h(z)) O(h(z)) O(1) O(1) !, with hdefined in (2.4). It will turn out that (y1)–(y4) is related to the inverse of the solution of (Y1)–(Y4). For any integrable N×Nmatrix-valued function Fon [a, b], C(F)(z) := 1 2πi Zb a F(t) t−zdt defines the Cauchy or Stieltjes transform of F, which is a matrix-valued and analytic function in C\[a, b]. Let us introduce the matrix polynomials of the second kind (of degree n−1), defined by Qn(x) = Zb a Pn(t)−Pn(x) t−xW(t)dt, n ≥0.(2.6) 6 F.A. Gr¨unbaum, M.D. de la Iglesia and A. Mart´ınez-Finkelshtein We have 2πiC(PnW)(z) = Qn(z)+2πiPn(z)C(W)(z), where 2πiC(W) is the m-function of the weight matrix W. Finally, if κnis the leading coefficient of any corresponding normalized polynomial Pn, we denote γn=κ∗ nκn.(2.7) Observe that a priori γndepends on the selection of κn. The following is a complete analogue of the well-known theorem from [28], although its proof requires some additional considerations: Theorem 2.2. The unique solution of the RHP (Y1)–(Y4) is Yn(z) = Rn(z)Y0(z), n ≥0,(2.8) where Y0(z) = INC(W)(z) 0IN, and the transfer matrix Rnis a matrix polynomial given by R0(z) = I, Rn(z) = κ−1 nPn(z) (2πiκn)−1Qn(z) −2πiκ∗ n−1Pn−1(z)−κ∗ n−1Qn−1(z)!, n ∈N.(2.9) Analogously, the unique solution of the RHP (y1)–(y4) is yn(z) = y0(z)rn(z), n ≥0,(2.10) where y0(z) = IN−C(W)(z) 0IN, and the transfer matrix rnis a matrix polynomial given by r0(z) = I, rn(z) = −Q∗ n−1(z)κn−1−(2πi)−1Q∗ n(z)κ−∗ n 2πiP∗ n−1(z)κn−1P∗ n(z)κ−∗ n!, n ∈N.(2.11) Moreover, for all n≥0, yn(z) = 0−IN IN0(Yn(z))T0IN −IN0= (Yn(z))−1,(2.12) det Yn(z)=1,for all z∈C.(2.13) Here ATdenotes the transpose of the matrix A. Remark 2.3. In the scalar case, when Ynis a 2 ×2 matrix, there is no need to consider simultaneously both RHPs (Y1)–(Y4) and (y1)–(y4) since the existence of (Yn)−1and (2.11), (2.12) follow directly from (2.13). In the general (2N)×(2N) case more care should be put in the analysis of the local behavior at the endpoints a,b, and (2.12) is no longer straightforward. Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization 7 Remark 2.4. The solution Ynof (Y1)–(Y4) satisfies the following symmetry relation: Yn(z) = IN0 0−INYn(¯z)IN0 0−IN,(2.14) which yields some obvious consequences for the block entries of Yn. This relation is established using the invariance of (Y1)–(Y4) by such a conjugation. Remark 2.5. Alternatively, we can write the solution (2.8) as Yn(z) = b Pn(z)C(b PnW)(z) −2πiγn−1b Pn−1(z)−2πiγn−1C(b Pn−1W)(z)!, n ∈N,(2.15) where b Pndenote the monic MOPRL of degree nfor the weight matrix W. In the same vein, (Yn)−1(z) =  −2πiC(Wb P∗ n−1)(z)γn−1−C(Wb P∗ n)(z) 2πi b P∗ n−1(z)γn−1b P∗ n(z) , n ∈N.(2.16) Remark 2.6. Although κnis defined up to a left unitary factor, the matrix coefficient γn in (2.7) is unique. This is a consequence of the uniqueness of the solutions of the RHP above. Remark 2.7. Formulas (2.8) and (2.10) show that generically the behavior of Ynand (Yn)−1 at the endpoints (and in general, any singular point) of the support of Wis given by the local behavior of the Cauchy transform of the orthogonality weight, also known as its m-function, see [50,§1.2]. Proof of Theorem 2.2.The fact that (2.8) or (2.15) is a solution of the RHP (Y1)–(Y4) is established following the proof for the scalar case, see e.g. [12], and taking into account that the orthogonality of b Pnin (1.1) is equivalent to the homogeneous system Zb a xjb Pn(x)W(x)dx =0, j = 0,1, . . . , n −1. The same applies to (2.10) or (2.16) and the RHP (y1)–(y4). Consider the function Zn(z) = Yn(z)yn(z); from (Y2) and (y2) it follows that it has no jump across (a, b), and by (Y4), (y4), Zn(z) = O(h(z)), z →c∈ {a, b}. Hence, Znhas only removable singularities at the finite endpoints of the support of the weight, and thus is an entire function. It remains to observe that Zn(∞) = I2Nto conclude that Zn(z) = I2Nfor all z∈C, which proves that yn(z) = (Yn(z))−1, as well as the uniqueness of both solutions. The first identity in (2.12) can be established by direct calculation or by observing that this transformation carries the RHP (Y1)–(Y4) to the RHP (y1)–(y4). Finally, the scalar function det Yn(z) is analytic across [a, b], and by (2.8), (2.9), det Yn(z) = det Rn(z) can have only removable singularities at the (finite) endpoints of the support of W. Hence, det Yn(z) is an entire function; since by (2.3), det Yn(∞) = 1, we conclude that it is identically 1.  8 F.A. Gr¨unbaum, M.D. de la Iglesia and A. Mart´ınez-Finkelshtein Proposition 2.8. Let Pn=κnb Pn,n≥0, be a sequence of orthonormal MOPRL, and let (Qn)n be the corresponding matrix polynomials of the second kind (2.6). For these polynomials, define An=κn−1κ−1 n.(2.17) Then Qn(z)P∗ n−1(z)−Pn(z)Q∗ n−1(z) = A−1 n,(2.18) Qn(z)P∗ n(z) = Pn(z)Q∗ n(z).(2.19) Moreover, P∗ n−1(x)AnPn(x)and Q∗ n−1(x)AnQn(x)are Hermitian for all x∈Rand n≥0. Remark 2.9. Identities (2.18) and (2.19) are also known as the Liouville–Ostrogradski formula and the Hermitian property, respectively. Both were originally derived for MOPRL in [20], using a different approach. Proof. Using the explicit expressions for Ynand (Yn)−1(written in the form (2.8) and (2.10)), from block entries (1,1) and (2,2) of the identity Yn(Yn)−1=Iwe obtain the Liouville– Ostrogradski formula (2.18) while from block entries (1,2) and (2,1) we get the so-called Hermitian property (2.19). Block entries (1,2) and (2,1) of (Yn)−1Yn=Iyield the commutativity relations P∗ n−1(z)AnPn(z) = P∗ n(z)A∗ nPn−1(z),Q∗ n−1(z)AnQn(z) = Q∗ n(z)A∗ nQn−1(z), which show that P∗ n−1(x)AnPn(x) and Q∗ n−1(x)AnQn(x) are Hermitian for all x∈Rand n≥0.  For the formulation of the following results we need to introduce a new set of parameters, bn,k =Zb a xkb Pn(x)W(x)dx, k =n, n + 1,..., (2.20) where b Pnare the monic MOPRL. In the spirit of [48], we can express them in terms of the coefficients an,j of the polynomials b Pn(see (1.1)) in a form suitable for numerical implementation: Lemma 2.10. Let Ωbe the block lower triangular matrix built up from the coefficients of the MOPRL (b Pn)n, Ω=     I a1,0I . . .. . .... an,0an,1· · · I     ,so that     b P0(x) b P1(x) . . . b Pn(x)     =Ω     I xI . . . xnI     ,I=IN. Then bn,k defined in (2.20)can be obtained from the last N×Nblock row of Ω−1as follows: b∗ n−k,nγn−k=Ω−1n+1,n−k+1, k = 0, . . . , n. (2.21) Proof. Observe that b∗ 0,n,...,b∗ n−1,n, b∗ n,n=Zb a xnWb P∗ 0,b P∗ 1,..., b P∗ ndx =Zb a xnW(I, xI, . . . , xnI)Ω∗dx = (µn,µn+1,...,µ2n)Ω∗, Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization 9 where µn=Rb axnW(x)dx are the moments of the weight matrix W. Analogously, by (1.2), (diag (γk)n k=0)−1= diag γ−1 kn k=0 =Zb a     b P0 b P1 . . . b Pn     Wb P∗ 0,b P∗ 1,..., b P∗ ndx =ΩZb a      I xI . . . xnI     WI, xI, . . . , xnIΩ∗=Ω     µ0µ1· · · µn µ1µ2· · · µn+1 . . .. . ..... . . µnµn+1 · · · µ2n     Ω∗. Hence, b∗ 0,nγ0,...,b∗ n−1,nγn−1,b∗ n,nγn =µn,µn+1,...,µ2n     µ0µ1· · · µn µ1µ2· · · µn+1 . . .. . ..... . . µnµn+1 · · · µ2n      −1 Ω−1=0,0,...,IΩ−1, which yields (2.21).  For what follows it is useful to single out the explicit expressions for bn−3,n,...,bn,n, that are obtained from Lemma 2.10 by direct computations: Corollary 2.11. Coefficients bn−k,n, for k= 0,1,2,3, are given by γnbn,n =I, γn−1bn−1,n =−a∗ n,n−1, γn−2bn−2,n =a∗ n−1,n−2a∗ n,n−1−a∗ n,n−2, γn−3bn−3,n =−a∗ n−2,n−3a∗ n−1,n−2a∗ n,n−1+a∗ n−1,n−3a∗ n,n−1+a∗ n−2,n−3a∗ n,n−2−a∗ n,n−3. Now we return to Theorem 2.2; as its immediate consequence we can relate the coefficients in the asymptotic expansion of Ynand of (Yn)−1, with the coefficients an,j and bn,j: Corollary 2.12. The coefficients Yn iin (2.3)are given by Yn i= an,n−i−1 2πibn,n+i−1 −2πiγn−1an−1,n−iγn−1bn−1,n+i−1 , i ≥0.(2.22) Analogously, the coefficients e Yn iin (2.5)are given by e Yn i=  bT n−1,n+i−1γn−1 1 2πibT n,n+i−1 2πiaT n−1,n−iγn−1aT n,n−i   , i ≥0.(2.23) Remark 2.13. One of the consequences of (2.14) is that an,j and bn,j appearing in (2.22), (2.23) belong to RN×N. The explicit expressions for coefficients of the asymptotic expansion above, combined with the obvious fact that Yn(Yn)−1=I2N, yield in a straightforward way the following identities: 16 F.A. Gr¨unbaum, M.D. de la Iglesia and A. Mart´ınez-Finkelshtein and for m∈N(see (2.3)), Xn(z) = I2N+ m X i=1 Yn i zi+O(1/zm+1)! znIN0 0z−nIN!V(z), z → ∞.(3.7) Since Xn(z) is invertible in C\R, we can consider the matrix function Fn(z)=d dz Xn(z)Xn(z)−1. Again, Fn(z) is analytic in C\Rand on the real line, (Fn)+(x)=(Fn)−(x), which implies that it is an entire matrix function. From (2.3) and (3.7), for z→ ∞, d dz Yn(z)[Yn(z)]−1=O1 z, and combining it with (2.3), (2.5), (3.5) and (3.7), we get for z→ ∞, Fn(z) = d dz Xn(z)Xn(z)−1 = I2N+ m X i=1 Yn i zi+O1/zm+1! G(z)0 0−G∗(z)! × I2N+ m X i=1 e Yn i zi+O1/zm+1!+O(1/z).(3.8) By Liouville’s theorem, the right hand side in (3.8) will coincide with its polynomial part in the expansion at infinity, and its degree is no greater than the degree of G. To be more precise, if we assume that the degree of Gis m∈Nand denote G(z) = m X j=0 Mjzj,f Mj= Mj0 0−M∗ j!, then after dropping the negative powers of zin Fn(z) = I2N+ m X i=1 Yn i zi+O1/zm+1! ×  m X j=0 f Mjzj  I2N+ m X k=1 e Yn k zk+O1/zm+1!+O(1/z) we obtain that Fn(z) = m X k=0   m X j=k j−k X i=0 Yn if Mje Yn j−i−k zk,with Yn 0=e Yn 0=I2N.(3.9) For instance, if m= 0, Fn(z) = f M0,(3.10) for m= 1, Fn(z) = f M1z+f M0+Yn 1f M1+f M1e Yn 1,(3.11) and for m= 2, Fn(z) = f M2z2+ ( f M1+Yn 1f M2+f M2e Yn 1)z+f M0+Yn 1f M1+f M1e Yn 1 +Yn 2f M2+Yn 1f M2e Yn 1+f M2e Yn 2.(3.12) Finally, using the explicit expressions (2.22) and (2.23) in (3.9), we arrive at the following Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization 17 Theorem 3.1. Under assumptions (3.2)and (3.4), with G(z) = m P j=0 Mjzj∈Pm, the matrix function Xndefined in (3.6)satisfies the following first-order differential equation with polynomial coefficients: d dz Xn(z) = Fn(z;G)Xn(z),(3.13) where Fn(z;G) =  −Bn(z;G)−1 2πiγ−1 nAn(z;G) 2πi An−1(z;G)γn−1B∗ n(z;G) ,(3.14) Anand Bnare matrix polynomials, An(z;G) = −γn  m−1 X j=0 bn,n+m−j−1∆∗ j,n(z) + ∆j,n(z)b∗ n,n+m−j−1 ,(3.15) Bn(z;G) = −  m X j=0 ∆j,n(z)b∗ n−1,n+m−j−1+bn,n+m−j−2∆∗ j,n−1(z) γn−1,(3.16) and the coefficients ∆j,n(z)are given by ∆j,n(z) = j X k=0 b Pn,k(z)Mm−j+k,b Pn,k(z) = zkI+an,n−1zk−1+· · · +an,n−k. Moreover, Fn(·;G)is linear in G: Fn(z;G1+G2) = Fn(z;G1) + Fn(z;G2).(3.17) Remark 3.2. The coefficients bn,k were introduced in (2.20) and discussed in Lemma 2.10. It is worth observing the similarity of this result with Theorem 2.17; in the present situation we can compute the differential equation for Xndirectly in terms of the coefficients of the MOPRL b Pn without considering integrals or studying the behavior at the endpoints. Finally, we have once again that γnA∗ n(z;G) = An(z;G)γn.(3.18) Again, for lowest degrees min Theorem 3.1 we can use (3.10), (3.11), (3.12), Lemma 2.10 and Corollary 2.11 in order to write the coefficient matrix (3.14) of the differential equation explicitly. For instance, for m= 0, Fn(z;G) = M00 0−M∗ 0; for m= 1, Fn(z;G) = G(z) + an,n−1M1−M1an,n−1 1 2πi(γ−1 nM∗ 1+M1γ−1 n) −2πi(γn−1M1+M∗ 1γn−1)−G∗(z) + a∗ n,n−1M∗ 1−M∗ 1a∗ n,n−1 ,(3.19) and for m= 2, with the notation (3.14), Bn(z;G) = −G(z)−(an,n−1M2−M2an,n−1)z+M1an,n−1−an,n−1M1−an,n−2M2 18 F.A. Gr¨unbaum, M.D. de la Iglesia and A. Mart´ınez-Finkelshtein −M2(an+1,nan,n−1+an+1,n−1) + an,n−1M2an,n−1−γ−1 nM∗ 2γn−1, An(z;G) = −M∗ 2z−M∗ 1+a∗ n+1,nM∗ 2−M∗ 2a∗ n,n−1 −γn(M2z+M1−M2an+1,n +an,n−1M2)γ−1 n. Observe from Theorem 3.1 that Bn(z;G) is a matrix polynomial of degree at most m and An(z;G) is a matrix polynomial of degree at most m−1. Finally, it is important to emphasize that in many situations we can exploit the implications of the freedom in the factorization (3.3) on the differential equation (3.13). This freedom, as we will see in Proposition 3.9, appears only in the matrix setting. Proposition 3.3. If under assumptions (3.2)and (3.4), with Ga polynomial, there exists a non-trivial matrix-valued function S, non-singular on C, smooth and unitary on R, such that H(z) = T(z)S0(z)S∗(z)T−1(z) (3.20) is also a polynomial, then e T=T S satisfies W(x) = e T(x)e T∗(x), x ∈R,e T0(z) = e G(z)e T(z), z ∈C, with e G(z) = G(z) + H(z). Moreover, the matrix Xndefined in (3.6), satisfies along with (3.13)–(3.16)the following relation: d dz Xn(z) = (Fn(z;G) + Fn(z;H)) Xn(z)−Xn(z) χ(z)0 0−χ∗(z)!,(3.21) with χ(z) = S0(z)S∗(z). Remark 3.4. Equations (3.13) and (3.21) are not necessarily trivially related, and in principle we could combine them in order to establish new relations for Xn, and thus, for Yn. Proof. The key observation is the formula (3.17), so that the coefficients Anand Bnin (3.13) depend linearly on G: An(·;G+H) = An(·;G) + An(·;H),Bn(·;G+H) = Bn(·;G) + Bn(·;H).(3.22)  3.2 Differential properties Taking into account the similarities between Theorems 2.17 and 3.1, we get immediately the analogue of Proposition 2.20 (the compatibility conditions): Proposition 3.5. Under assumptions (3.2)and (3.4), with G(z)a polynomial, the coefficients of the differential equation (3.13)satisfy the following recurrence relations: for every n≥0, I+Bn+1(z;G)(zI−αn)−(zI−αn)Bn(z;G) = A∗ n+1(z;G)βn+1 −βnA∗ n−1(z;G) (3.23) and Bn+1(z;G) + γ−1 nB∗ n(z;G)γn= (zI−αn)A∗ n(z;G),(3.24) where αnand βnare the coefficients of the three-term recurrence relation (2.30). As before, the ladder operators (the most basic differential properties for MOPRL) can be easily obtained by analyzing the first block column of Xnin the differential equation (3.13): Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization 19 Proposition 3.6. Under assumptions (3.2)and (3.4), with G(z)a matrix polynomial, the monic MOPRL (b Pn)nsatisfy the following difference-differential relations (lowering and raising operators, respectively): b P0 n(z) + b Pn(z)G(z) = −Bn(z;G)b Pn(z) + A∗ n(z;G)βnb Pn−1(z) (3.25) and b P0 n(z) + b Pn(z)G(z) = A∗ n(z;G)(zI−αn)− Bn(z;G)b Pn(z)− A∗ n(z;G)b Pn+1(z),(3.26) where Anand Bnare given by (3.15)and (3.16), respectively. Proof. Block entry (1,1) of (3.13) gives the lowering operator (3.25), using (2.31) and (3.18), while block entry (2,1) of (3.13) gives the raising operator (3.26) using (3.24) and (3.18).  Remark 3.7. Comparing with the results of Corollary 2.21, notice that these ladder operators contain a term with the MOPRL multiplied on the left. In the situation described in Proposition 3.3, we can in principle exploit the non-uniqueness of the relations above in order to derive further relations for the MOPRL. The following result shows the possibility of the existence of 2-terms recurrence relations for the family b Pn, a phenomenon that has not been reported before in the theory of MOPRL. Corollary 3.8. Under conditions of Proposition 3.3, the family of monic MOPRL (b Pn)nsatisfies, along with (3.25)–(3.26), the following relations: b Pn(z)H(z) = −Bn(z;H)b Pn(z) + A∗ n(z;H)βnb Pn−1(z) (3.27) and b Pn(z)H(z) = A∗ n(z;H)(zI−αn)− Bn(z;H)b Pn(z)− A∗ n(z;H)b Pn+1(z),(3.28) where His defined in (3.20). In particular, provided An(z;G)is invertible for all z∈C, b Pn(z)H(z) + Bn(z;H)b Pn(z)− A∗ n(z;H)A−∗ n(z;G) ×b P0 n(z) + b Pn(z)G(z) + Bn(z;G)b Pn(z)=0.(3.29) Proof. To get (3.27), we subtract (3.13) and (3.21) and then evaluate the (1,1) block entry. (3.28) is obtained from (3.27) and the three term recurrence relation. Finally, (3.29) is obtained by replacing (3.25) in (3.27).  Equations (3.27) and (3.28) are known as the 0-th order ladder operators, while (3.29) is a first order differential relation for the MOPRL. In some situations they yield trivial identities; this is always true in the scalar case, as the following proposition shows: Proposition 3.9. Assume that under the conditions of Proposition 3.3,χ(z) = S0(z)S∗(z) = ip(z)I, where pis a scalar polynomial of degree m. Then An(z;H) = 0,and Bn(z;H) = −ip(z)I. Proof. From (3.20), H=T χT −1, and by linearity of Anand Bnin H(see (3.22)), it is enough to prove the formulas above for monomials zkI. In this case, formulas (3.15) and (3.16) simplify considerably. Using Proposition 2.14 gives An(z;H) = 0and Bn(z;H) = −izkI. 20 F.A. Gr¨unbaum, M.D. de la Iglesia and A. Mart´ınez-Finkelshtein In the scalar case (N= 1) the only smooth and unitary function s(x) on Rhas the form s(z) = eip(x), with preal-valued, and s0(z)s(¯z) = ip(z), so that the assumption that pis a polynomial brings us to the situation described in Proposition 3.9. This explains why in the scalar setting we never get nontrivial 0-th order ladder operators. In the scalar case one cannot have first order differential relations for the MOPRL, such as those given in (3.29). The general matrix case is much richer and complex, as will be illustrated with some examples in the next section. We finish this section with a class of second-order differential equations satisfied by the MOPRL ( b Pn)n. As a consequence of the Proposition 3.6 we have the following Proposition 3.10. Under assumptions (3.2)and (3.4), with G(z)a polynomial, the MOPRL (b Pn)nsatisfy the following second-order differential equation b P00 n+ 2 b P0 nG+b PnG0+G2+Mnb P0 n+Nnb Pn+Mnb PnG=0,(3.30) where An(z) = An(z;G),Bn(z) = Bn(z;G), Mn(z) = Mn(z;G) =−(A∗ n(z))0A−∗ n(z) + Bn(z)− A∗ n(z)(zI−αn) + A∗ n(z)Bn+1(z)A−∗ n(z), and Nn(z) = Nn(z;G) = Mn(z)Bn(z)− B2 n(z) + B0 n(z) + A∗ n(z)βnA∗ n−1(z), provided that the inverse of An(z;G)exists for z∈C. Here, again, −∗ denotes the conjugate transpose of the inverse. Proof. Differentiate the lowering operator (3.25) and substitute the raising operator (3.26) evaluated at n−1.  Remark 3.11. Note that we can easily obtain another second-order differential operator satisfied by the orthogonal polynomials reversing the ladder operators. This new differential equation needs not be in principle the same as (3.30). Nevertheless, it is straightforward to see that both equations are equivalent using the compatibility conditions (3.23) and (3.24). The equation (3.30) does not have the form of the right hand side differential operator considered for instance in [23], due to the terms Mnb P0 n,Nnb Pnand Mnb PnG. In some cases, under additional assumptions on the weight, (3.30) can be reduced further, as we will see in the following section. 4 Illustrative examples In this section we study a number of examples of weights Wwhich are smooth and non-vanishing on the whole real line; this assumption simplifies the Riemann–Hilbert formulation because in this case one does not consider the local conditions (Y4) (see Section 2.1). Additionally, without loss of generality, we take W(0) = I. For convenience, we consider the weights of the form W(x) = e−2q(x)U(x)U∗(x), x ∈R,(4.1) Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization 21 where qis a scalar real-valued function, so that with the notation (3.2) and (3.4), we may take T(x) = e−q(x)U(x) and G(x) = −q0(x)I+U0(x)U−1(x).(4.2) We are interested in the case when Gis a matrix polynomial. Hence, keeping up with previous hypotheses, we will assume that qis a (scalar) polynomial of even degree with real coefficients and a positive leading coefficient, and that U0(x)U−1(x) is a matrix polynomial. Since our main goal here is to generate a set of examples (some new, some already well known), we restrict the degree of qto either 2 (the Hermite case) or 4 (the Freud case), and the degree of U0(x)U−1(x) to at most 1. We start by considering the case when U0(x)U−1(x) is a monomial, i.e. either U0(x)U−1(x) = Aor U0(x)U−1(x)=2Bxfor constant matrices A,B∈CN×N. Taking into account the linearity (3.22) of the coefficients Anand Bnin (3.15), (3.16), we can obtain the differential equation (3.13) for the general case of U0(x)U−1(x) = A+ 2Bx. However, finding the corresponding orthogonality weight is more involved: when Aand Bdo not commute, solving U0(x)U−1(x) = A+ 2Bxis not straightforward. In Section 4.1.3 we will discuss some examples, which yield explicit expressions for U(x), related to this case. Recently an example has been found in [3] where a weight matrix supported in the real line is explicitly given (but not of the type (4.1)), when G(x) is a matrix polynomial of degree N with in general non-commuting coefficients. 4.1 The Hermite case For q(x) = x2/2, let us consider two cases, first U0(x)U−1(x) = Aand then U0(x)U−1(x) = 2Bx. We end up this Section by discussing briefly the case of U0(x)U−1(x) = A+ 2Bx. 4.1.1 U0(x)U−1(x) = A The differential equation U0(x)U−1(x) = A, so that U(x) = eAx,T(x) = e−x2/2eAx, and the weight matrix (4.1) is given by W(x) = e−x2eAxeA∗x,A∈CN×N, x ∈R.(4.3) The matrix Ghas the form G(x) = −xI+A, and according to (3.19), An(x;G) = 2I,Bn(x;G) = −G(x). The compatibility conditions of Proposition 3.5 yield 2(βn+1 −βn) = I+Aαn−αnA,(4.4) and αn=1 2(A+γ−1 nA∗γn).(4.5) The lowering and raising operators from Proposition 3.6 are reduced now to b P0 n(x) + b Pn(x)A−Ab Pn(x) = 2βnb Pn−1(x),(4.6) and −b P0 n(x)+2xb Pn(x) + Ab Pn(x)−b Pn(x)A−2αnb Pn(x) = 2 b Pn+1(x). 22 F.A. Gr¨unbaum, M.D. de la Iglesia and A. Mart´ınez-Finkelshtein Summing up the telescopic relation in (4.4) (or comparing the O(xn−1) term in (4.6)) gives βn=1 2(nI+an,n−1A−Aan,n−1).(4.7) With the notation of Proposition 3.10, Mn(x) = 2 (αn−A),Nn(x) = −2 (αn−A)G(x)−G2(x) + I+ 4βn, so that the differential equation (3.30) for the monic polynomials b Pn, orthogonal with respect to the weight (4.3), boils down to b P00 n(x)+2b P0 n(x)(A−xI) + b Pn(x)(A2−2xA) = (−2xA+A2−4βn)b Pn(x) + 2(A−αn)( b P0 n(x) + b Pn(x)A−Ab Pn(x)).(4.8) These formulas hold for any constant matrix A, and we cannot expect important simplifications without narrowing the class of the weights further. This can be done assuming in addition that the hypotheses of Proposition 3.3 hold. This, as it was shown in [23,24], imposes additional constraints on the weight W. Since the construction is described in detail in [23,24], the exposition in this part will be rather sketchy. The matrix Hfrom (3.20), restricted to R, can be written as H(x) = T(x)χT −1(x) = eAxχe−Ax =χ+ adA(χ)x+ ad2 A(χ)x2 2+· · · =X k≥0 adk A(χ)xk k!.(4.9) Here χ(x) = S0(x)S∗(x) is skew-Hermitian on R, adAis the commutator given by adA(χ) = Aχ −χA, and we define recursively ad0 A(χ) = χ,adn+1 A(χ) = adA(adn A(χ)) for n≥1. The simplest situation obtains when the right hand side in (4.9) is constant; as it was shown in Lemma 2.4 of [24], this assumption yields χ=iaIfor certain a∈R(which is the case discussed in Proposition 3.9), when there are no new ladder operators. Consequently, for the first non-trivial situation we must assume that (4.9) is a matrix polynomial of degree at least one. We consider two situations that yield degree exactly one (see [23,24] for motivations and further details). Let us define a nilpotent matrix of the form L= N−1 X k=1 νkEk,k+1, νk∈C\ {0},(4.10) where Eij is a matrix with 1 at entry (i, j) and 0 elsewhere, and a diagonal matrix J= N X k=1 (N−k)Ek,k.(4.11) For the first non-trivial example we assume that A=Land χ=iJ, so that adA(χ) = −A(4.12) and ad2 A(χ) = 0. Thus, H(x) = eAxχe−Ax=i(J−Ax), and by (3.19), An(x;H) = i(−A∗+γnAγ−1 n)=2i(α∗ n−A∗), Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization 23 Bn(x;H) = i(Ax−J+an,n−1A−Aan,n−1) = i(Ax−J+ 2βn−nI). From the compatibility conditions of Proposition 3.5 we get Jαn−αnJ+αn=A+1 2A2αn−αnA2and J−γ−1 nJγn=Aαn+αnA−2α2 n,(4.13) where we have used relations (4.4), (4.5), (4.7) and (A−2αn)βn=βn(A−2αn−1), which can be easily proved using (2.31) and (4.5). Thus, from Corollary 3.8, (4.4) and (4.13), the lowering and raising operators (of the 0-th order) are b Pn(x)J−Jb Pn(x)−x(b Pn(x)A−Ab Pn(x)) + 2βnb Pn(x)−nb Pn(x) = 2(A−αn)βnb Pn−1(x), and b Pn(x)(J−xA)−γ−1 n(J−xA∗)γnb Pn(x)+2βn+1 b Pn(x)−(n+ 1) b Pn(x) = 2(αn−A)b Pn+1(x), respectively, which yields by (3.29) the first order relation (A−αn)b P0 n(x)+(A−αn+xI)( b Pn(x)A−Ab Pn(x)) −2βnb Pn(x) =b Pn(x)J−Jb Pn(x)−nb Pn(x).(4.14) These identities hold in addition to the second order equation obtained in (4.8), valid as we recall, for any A. As far as we are aware of, these relations are new. We use them to simplify the differential equation remarkably: multiplying (4.14) by 2 and plugging it into (4.8) gives us for A=L, b P00 n(x)+2b P0 n(x)(A−xI) + b Pn(x)A2−2J=−2nI+A2−2Jb Pn(x), which is a linear second-order differential equation with coefficients in the right hand side independent on n(a.k.a. Sturm–Liouville equation with polynomial coefficients), studied by Dur´an and Gr¨unbaum in [23]. In the second case we take A=L(I+L)−1= N−1 X j=1 (−1)j−1Lj,χ=iJ, where Lis given in (4.10) and Jin (4.11). Consequently, adA(χ) = −A+A2,ad2 A(χ) = 0,(4.15) and H(x) = eAxχe−Ax=iJ−A−A2x. In this case, An(x;H) = i−A∗+ (A∗)2+γnA−A2γ−1 n 24 F.A. Gr¨unbaum, M.D. de la Iglesia and A. Mart´ınez-Finkelshtein = 2i(α∗ n−A∗−(α∗ n−A∗)α∗ n−α∗ n(α∗ n−A∗)) and Bn(x;H) = iA−A2x−J+an,n−1A−A2−A−A2an,n−1 =iA−A2x−J+ 2βn−nI−2(Aβn+βnA)+2nA. As to new results, we get the lowering operator, b Pn(x)J−xA−A2−J−xA−A2b Pn(x) = (−2βn+n(I−2A) + 2(Aβn+βnA)) b Pn(x) −2(αn−A−(αn−A)αn−αn(αn−A))βnb Pn−1(x), the raising operator b Pn(x)J−xA−A2−J−xA−A2b Pn(x) = (−2βn+n(I−2A) + 2(Aβn+βnA)) b Pn(x) −2(αn−A−(αn−A)αn−αn(αn−A))((xI−αn)b Pn(x)−b Pn+1(x)), and the first-order differential equation (3.29): (A−αn−(A−αn)αn−αn(A−αn)) b P0 n(x) =b Pn(x)J−Jb Pn(x)−xb Pn(x)A2−A2b Pn(x) + (xI−A+αn+ (A−αn)αn+αn(A−αn))( b Pn(x)A−Ab Pn(x)) + (2βn+n(2A−I)−2(Aβn+βnA)) b Pn(x). Multiplying this equation by 2 and plugging it into the second-order differential equation (4.8) gives b P00 n(x)+2b P0 n(x)(A−xI) + b Pn(x)A2−2xA2−2J=A2−2xA2−2Jb Pn(x) + (2n(2A−I)−4(Aβn+βnA) + 2((αn−A)αn+αn(αn−A))A)b Pn(x) −2((αn−A)αn+αn(αn−A))( b P0 n(x) + b Pn(x)A). This differential equation is not of Sturm–Liouville type considered by Dur´an and Gr¨unbaum in [23], but nevertheless we have been able to give a number of differential equations of first and second order satisfied by MOPRL with respect to a weight matrix that have not been considered up to this point. 4.1.2 U0(x)U−1(x)=2Bx The weight matrix (4.1) is now given by W(x) = e−x2eBx2eB∗x2,B∈CN×N, x ∈R, where we assume that Bis chosen such that all the moments exist. By (4.2), T(x) = e−x2/2eBx2and G(x) = (2B−I)x. The weight matrix is an even matrix function, so that all moments of odd order vanish. As a consequence, an,n−1=0and αn=0. Then An(x;G) = 2I−B∗−γnBγ−1 n,Bn(x;G)=(I−2B)x, Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization 25 and there will be only one compatibility condition: 2I−B−γ−1 n+1B∗γn+1βn+1 −2βnI−B−γ−1 n−1B∗γn−1=I. From (4.16) we get the explicit expression for βnvia 2I−B−γ−1 nB∗γnβn=nI+ 2(an,n−2B−Ban,n−2). The lowering and raising operators are b P0 n(x)+2x(b Pn(x)B−Bb Pn(x)) = 2I−B−γ−1 nB∗γnβnb Pn−1(x),(4.16) and b P0 n(x)+2xb Pn(x)B−Bb Pn(x)= 2I−B−γ−1 nB∗γnxb Pn(x)−b Pn+1(x), respectively. This yields the following second-order differential equation b P00 n(x)+2xb P0 n(x)(2B−I)+2xγ−1 nB∗γn−Lnb P0 n(x)+4x2b Pn(x)B2−B2b Pn(x) + 4Knb Pn(x) + 2−4x2I+ 4x2γ−1 nB∗γn−Lnb Pn(x)B−Bb Pn(x)=0, where Ln=I−B−γ−1 nB∗γnBI−B−γ−1 nB∗γn−1, and Kn=I−B−γ−1 nB∗γnβnI−B−γ−1 n−1B∗γn−1. These expressions are valid for any B(as long as all moments of the weight matrix exist). For further simplifications we can assume again that the hypotheses of Proposition 3.3 hold, and consider the cases given by the algebraic relations (4.12) and (4.15). For instance, when (4.15) holds, we obtain again a second-order differential equation with coefficients in the right hand side independent on n, studied by Dur´an and Gr¨unbaum in [23]. The computations are similar and will be omitted for the sake of brevity. 4.1.3 Other cases Assume now that we have a linear combination of the previous two cases, i.e. U0(x)U−1(x) = A+ 2Bx, A,B∈CN×N,with U(0) = I,(4.17) in which case the matrix G(x) is given by G(x) = A+ (2B−I)x. This is all that is needed to calculate the coefficients An(x;G) and Bn(x;G) and consequently the compatibility conditions, the ladder operators and the differential relations. However, for the corresponding orthogonality weight we need to solve (4.17) explicitly, which may be non-trivial (unless Aand Bcommute). In general, this solution can be given in terms of the time ordered exponential U(x) = : eRx 0(A+Bs)ds :,A,B∈CN×N, x ∈R. This is obtained by rewriting the differential equation as an integral one and “solving” it by iteration. In general, one cannot give an explicit expression for this infinite sum. An example of explicitly solvable non-trivial equation (4.17) can be found in [22, Theorem 1.1], where non-commuting Aand Bare constructed (as a linear combination of the matrices Land J