Symmetry, Integrability and Geometry: Methods and Applications Properties of Matrix Orthogonal Polynomials via their Riemann–Hilbert Characterization

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.


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 ≤ n in C N ×N and a scalar variable x can be defined as an expression of the form where A j 's are constant matrices in C N ×N . In what follows we consider N fixed, and denote by P n the family of matrix polynomials of degree ≤ n in C N ×N , as well as by P := n≥0 P n . We use preferably boldface letters to denote matrices, and standard font for scalars. We also use I N for the N × N identity 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 = (W ij ) : (a, b) → GL(N, R), defined and positive definite on a finite or infinite interval (a, b) ⊂ R. We will assume that all W ij and W ′ ij have finite moments: where the integration of a matrix function is applied entry-wise. For any two P , Q ∈ P, the weight W induces two matrix-valued "inner products", 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 P W := (Tr P , P W ) 1/2 , and assume that W is non-trivial, in the sense that P W > 0 for every non-zero matrix polynomial P . In this case (see [11,Lemma 2.3] or [51,Proposition 4.2.3]), (P , P ) W is 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 ( P n ) n of monic orthogonal polynomials such that P 0 = I, P n (x) = x n I + n−1 j=0 a n,j x j , ( P n , Q) W = 0 for every Q ∈ P n−1 , n ∈ N, (1.1) as well as matrix polynomials (P n ) n , "orthonormal" with respect to W , such that P n (x) = κ n P n , (P n , P m ) W = δ n,m I. (1.2) 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 3 we 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 4 we narrow the analysis of Sections 2 and 3 to 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, 3 and in [4,Section 2], the focus of both contributions is different, and in this sense, complementary.

.1 Formulation and basic properties
In this section we discuss the Riemann-Hilbert problem (RHP) related to MOPRL with respect to a N × N weight matrix W supported on an interval [a , b ] ⊂ R; this interval can be either bounded or unbounded. We assume W continuous and non-vanishing on (a, b), and that at any finite endpoint of the support the weight W has at worse a power-type singularity, that is, W is of the form where c ∈ {a, b}, c = ±∞, and W is 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 × 2N matrices 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 W consists in finding a matrix function (Y2) Y n has on (a, b) continuous boundary values Y n + (resp., Y n − ) from the upper (resp., lower) half plane, such that (where the asymptotic term O(1/z m+1 ) depends on n). Along with the RH problem (Y1)-(Y4) we can consider the following dual problem: finding a matrix function y n : (y3) As z → ∞, for every m ∈ N we have with h defined in (2.4).
It will turn out that (y1)-(y4) is related to the inverse of the solution of (Y1)-(Y4).
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 where 2πiC(W ) is the m-function of the weight matrix W . Finally, if κ n is the leading coefficient of any corresponding normalized polynomial P n , we denote Observe that a priori γ n depends 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: and the transfer matrix R n is a matrix polynomial given by R 0 (z) = I, Analogously, the unique solution of the RHP (y1)-(y4) is y n (z) = y 0 (z)r n (z), n ≥ 0, (2.10) where and the transfer matrix r n is a matrix polynomial given by r 0 (z) = I, Moreover, for all n ≥ 0, Here A T denotes the transpose of the matrix A.
Remark 2.3. In the scalar case, when Y n is a 2 × 2 matrix, there is no need to consider simultaneously both RHPs (Y1)-(Y4) and (y1)-(y4) since the existence of (Y n ) −1 and (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.
Remark 2.4. The solution Y n of (Y1)-(Y4) satisfies the following symmetry relation: which yields some obvious consequences for the block entries of Y n . 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 where P n denote the monic MOPRL of degree n for the weight matrix W . In the same vein, Remark 2.6. Although κ n is 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 Y n and (Y n ) −1 at the endpoints (and in general, any singular point) of the support of W is given by the local behavior of the Cauchy transform of the orthogonality weight, also known as its m-function, see 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 P n in (1.1) is equivalent to the homogeneous system b a x j P n (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 Z n (z) = Y n (z)y n (z); from (Y2) and (y2) it follows that it has no jump across (a, b), and by (Y4), (y4), Hence, Z n has only removable singularities at the finite endpoints of the support of the weight, and thus is an entire function. It remains to observe that Z n (∞) = I 2N to conclude that Z n (z) = I 2N for all z ∈ C, which proves that y n (z) = (Y n (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 Y n (z) is analytic across [a, b], and by (2.8), (2.9), det Y n (z) = det R n (z) can have only removable singularities at the (finite) endpoints of the support of W . Hence, det Y n (z) is an entire function; since by (2.3), det Y n (∞) = 1, we conclude that it is identically 1.
Proposition 2.8. Let P n = κ n P n , n ≥ 0, be a sequence of orthonormal MOPRL, and let (Q n ) n be the corresponding matrix polynomials of the second kind (2.6). For these polynomials, define Moreover, P * n−1 (x)A n P n (x) and Q * n−1 (x)A n Q n (x) are Hermitian for all x ∈ R and 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 Y n and (Y n ) −1 (written in the form (2.8) and (2.10)), from block entries (1, 1) and (2, 2) of the identity Y n (Y n ) −1 = I we 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 (Y n ) −1 Y n = I yield the commutativity relations which show that P * n−1 (x)A n P n (x) and Q * n−1 (x)A n Q n (x) are Hermitian for all x ∈ R and n ≥ 0.
For the formulation of the following results we need to introduce a new set of parameters, b n,k = b a x k P n (x)W (x) dx, k = n, n + 1, . . . , (2.20) where P n are the monic MOPRL. In the spirit of [48], we can express them in terms of the coefficients a n,j of the polynomials P n (see (1.1)) in a form suitable for numerical implementation: Let Ω be the block lower triangular matrix built up from the coefficients of the MOPRL ( P n ) n , . .
Now we return to Theorem 2.2; as its immediate consequence we can relate the coefficients in the asymptotic expansion of Y n and of (Y n ) −1 , with the coefficients a n,j and b n,j : Analogously, the coefficients Y n i in (2.5) are given by Remark 2.13. One of the consequences of (2.14) is that a n,j and b n,j appearing in (2.22), (2.23) belong to R N ×N .
The explicit expressions for coefficients of the asymptotic expansion above, combined with the obvious fact that Y n (Y n ) −1 = I 2N , yield in a straightforward way the following identities: Proposition 2.14. The coefficients of the monic MOPRL P n and the coefficients of the Cauchy transform C( P n W ) (i.e. the sequences b n,k ) satisfy the following relations m j=0 a n,n−m+j b * n−1, In particular, for m = 1 we recover the second identity of Corollary 2.11.
We finish this section with the matrix analogue of the Christoffel-Darboux (CD) formula for MOPRL, written in terms of the solution of the characterizing RHP. For an orthonormal family (P n ) n we consider the kernel Then the CD formula (see [20]) reads as Proposition 2.15. The matrix CD kernel (2.26) satisfies, for x, y ∈ R, The proof follows taking into account the Christoffel-Darboux formula (2.27), the expression for A n in (2.17) and the expression of the inverse (Y n ) −1 in (2.16). The last proposition for general weight matrices of rectangular size can be found in [18].

Difference and differential relations
Both the MOPRL and the solution of the RHP given in Theorem 2.2 depend on two variables, the discrete n and the continuous z. We can derive further properties by analyzing the variation of this solution with respect to either variable, which yields linear relations of the form The general methodology to get these equations may be traced back to the original work of Gelfand, Levitan and other authors. The idea is that if the jump matrix for a RHP is independent of a variable, then a variation with respect to that variable leads to an identity. The fact that the jump matrix in (2.2) is independent of n allows one to obtain immediately the first identity in (2.28), which in turn is connected to the well-known three-term recurrence relation for MOPRL.
Theorem 2.16. There is a unique sequence of matrix coefficients α n such that the solution of the RHP Y n (z) satisfies the following first-order difference equation Consequently, the monic MOPRL P n satisfy the following three-term recurrence relation: The coefficients α n and β n , as well as γ n (defined in (2.7)) can be expressed either in terms of the elements of the solution of the RHP or in terms of the coefficients a n,j (see (1.1)) as follows: The behavior at the exceptional points gives that the singularities at a and b are removable, and R is an entire function. From (2.3) and Liouville's theorem we conclude that This proves (2.29) along with the expression for α n in terms of the solution of the RHP. The three term recurrence relation (2.30) for the MOPRL is obtained by considering the (1, 1) block entry of (2.32). Finally, the rest of the expressions for the coefficients is found using Corollary 2.12 (formula (2.22)).
We turn now to the dependence of Y n on the continuous variable z; notice however that the jump (2.2) along the real line does depend on this variable, so the application of the paradigm explained at the beginning of this subsection is not straightforward. At this stage we can only aspire to find a differential relation like in (2.28), but with a matrix F n depending upon the entries of Y n (or equivalently, upon the MOPRL themselves).
Theorem 2.17. The solution of the RHP for Y n (z) satisfies the following first-order matrix differential equation

35)
provided that all the functions in the right-hand sides of (2.34) and (2.35) exist for z ∈ C \ [a, b].
In particular, Remark 2.18. Formulas (2.34) and (2.35) were introduced in the scalar case in [7] and in the matrix case (for a normalized family) in [25]. In general, they bear a formal character.
In particular, the boundary terms in (2.34) and (2.35) may be undefined, depending on the behavior of W at the endpoints (see (2.1)). As it will be shown in the next section, these formulas simplify considerably after considering a special transformation of the RHP for Y n . Finally observe that the coefficients (2.34) and (2.35) are unique since Y n is unique.
. From the explicit expressions (2.15) and (2.16) we have that the block entries of R n are In particular, −(R n (z)) 11 = (R n (z)) * 22 and −2πiγ n (R n (z)) 12 = 1 2πi (R n+1 (z)) 21 γ −1 n . We will now derive formulas (2.34) and (2.35). For that purpose we need the following technical observations: (ii) For any differentiable and integrable matrix function F bounded at z = a, b, Identities in (i) follow by adding and subtracting C(P W P * n )(z) and using the orthogonality of the MOPRL, and (ii) is obtained by integration by parts and using that d We return to the proof of Theorem 2.17, showing in detail how to obtain B n (z). For A n (z) the computations are similar and will be omitted for the sake of brevity.
Using (i) from the previous lemma and the product rule of differentiation we get Applying now (ii) to the second term and canceling we get which yields (2.35).
Let us discuss now two main consequences of the existence of the first order matrix-valued differential equation (2.33).
One one hand, equations (2.28), known as the Lax pair for the RHP, are clearly overdetermined, so compatibility conditions (via cross-differentiation of both equations) yield also known as string equations. They allow us, for instance, to recover the sequence {F n } from the sequence {E n } and the initial value F 1 . These ideas relate also the RHP to nonlinear problems such as the nonlinear Schrödinger equation (NLS), see [17], and the Toda flow, see [27]. In our situation, the compatibility conditions (2.37) considered entry-wise imply the following: Proposition 2.20. The recurrence relations hold for every n ≥ 0, where α n and β n are the coefficients of the three-term recurrence relation (2.30).
Another consequence of Theorem 2.17 is the existence of the ladder operators for MOPRL: Corollary 2.21. The monic MOPRL ( P n ) n satisfy the following difference-differential relations (lowering and raising operators, respectively): and  [25]; the RH approach gives an alternative proof of these results. The ladder operators are the basic differential relations for MOPRL. It is well-known that they can be combined to build a secondorder differential equation satisfied by the polynomials. This fact was mentioned in [25], but no explicit expression of the differential equation was given; here we present it for completeness: Corollary 2.22. The polynomials ( P n ) n satisfy the following second-order differential equation and provided that the inverse of A * n (z) exists for z ∈ C \ [a, b]. Here and below, the superscript − * denotes the conjugate transpose of the inverse.
Proof . Differentiate the lowering operator (2.40) and substitute the raising operator (2.41) evaluated at n − 1.
Remark 2.23. 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 in principle be the same as (2.42). Nevertheless, it is straightforward to see that both equations are equivalent using the compatibility conditions (2.38) and (2.39).
Although Theorem 2.17 gives explicit expressions of A n (z) and B n (z), these coefficients are usually difficult to calculate. They are not even defined on the interval [a, b]. In the next section we will introduce some additional assumptions that simplify the differential equation considerably.
3 Transformation of the RHP when supp(W ) = R As it was pointed out in the previous section, the independence of the jump in (2.2) with respect to the discrete variable n allows one to obtain the three-term recurrence relation in a straightforward way. Under additional assumptions on the jumps we can perform a transformation of the RHP in such a way that a new jump is independent of the continuous variable z. This will have consequences on the resulting differential relations.
In the matrix case there is some extra freedom absent in the scalar situation, that gives us a whole new family of differential relations, and consequently, a whole class of ladder operators, some of them of the 0-th order, something that is not possible in the scalar case. The combination of the ladder operators will give rise to a class of second-order differential equations, some of them of order less than 2.
For simplicity, we will focus here on the case when [a, b] = R, assuming along this section that W is smooth and positive definite on the whole real axis R. Some ideas about how to handle the case of finite endpoints of the support of W are briefly explained in the final Section 5.

The transformation
As we have seen in Section 2.2, the recurrence relation (first identity in (2.28)) is a general fact intrinsic to the orthogonality of the polynomials. Simple differential relations are much more demanding, and we must impose at this stage further conditions on the orthogonality weight W .
Our immediate goal is to obtain an invertible transformation Y n → X n such that X n has a constant jump across R. We will consider where V is a matrix-valued function, analytic in C \ R and continuous up to R, and invertible for all z ∈ C. Then the jump matrix for X n on R is Observe that this kind of transformations does not affect the first difference equation in the Lax pair (2.28), but constant jumps will allow us to use the strategy of variation of the problem along z. Going down the path of simplification, we can try a block-diagonal matrix

Properties of Matrix Orthogonal Polynomials via their Riemann-Hilbert Characterization 15
where T is an invertible N × N matrix function. Then (3.1) boils down to These preliminary considerations motivate to consider a factorization of the weight in the form where T is a smooth matrix-valued function on R. Under our assumptions the existence of such a T is guaranteed, but not its uniqueness: performing its RQ decomposition it can be written as where T (x) is an upper triangular matrix and S(x) is an arbitrary smooth and unitary matrix (for each x ∈ R). This representation can go further taking into account that any unitary matrix can be written as S(x) = e Q(x) , where Q(x) is an skew-Hermitian matrix function. Additionally, since any skew-Hermitian matrix is normal, it has a factorization Q( where U is again unitary and D a diagonal matrix with real entries. Therefore T from (3.3) can be written as We narrow the choice of T (and W ) by imposing the additional constraint that there exists a matrix polynomial G such that in a neighborhood of the origin, provides an analytic extension of W to the whole plane. Following our previous discussion, we define the entire matrix-valued function Observe that V is invertible for all z ∈ C. By differentiating the identity T −1 T = I and using (3.4) we conclude that If Y n (z) is the solution of the original RHP, let us define Then straightforward computations show that X n (z) is analytic in C \ R, satisfies the jump condition X n + (x) = X n − (x) and for m ∈ N (see (2.3)), Since X n (z) is invertible in C\R, we can consider the matrix function F n (z)= d dz X n (z) X n (z) −1 . Again, F n (z) is analytic in C \ R and on the real line, (F n ) + (x) = (F n ) − (x), which implies that it is an entire matrix function. From (2.3) and (3.7), for z → ∞, and combining it with (2.3), (2.5), (3.5) and (3.7), we get for z → ∞, 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 G is m ∈ N and denote then after dropping the negative powers of z in For instance, if m = 0, for m = 1, and for m = 2, Finally, using the explicit expressions (2.22) and (2.23) in (3.9), we arrive at the following M j z j ∈ P m , the matrix function X n defined in (3.6) satisfies the following first-order differential equation with polynomial coefficients:

14)
A n and B n are matrix polynomials, 16) and the coefficients ∆ j,n (z) are given by ∆ j,n (z) = j k=0 P n,k (z)M m−j+k , P n,k (z) = z k I + a n,n−1 z k−1 + · · · + a n,n−k .
Moreover, F n (·; G) is linear in G: Remark 3.2. The coefficients b n,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 X n directly in terms of the coefficients of the MOPRL P n without considering integrals or studying the behavior at the endpoints. Finally, we have once again that γ n A * n (z; G) = A n (z; G)γ n .  and for m = 2, with the notation (3.14), B n (z; G) = −G(z) − (a n,n−1 M 2 − M 2 a n,n−1 )z + M 1 a n,n−1 − a n,n−1 M 1 − a n,n−2 M 2 − M 2 (a n+1,n a n,n−1 + a n+1,n−1 ) + a n,n−1 M 2 a n,n−1 − γ −1 n M * 2 γ n−1 , A n (z; G) = −M * 2 z − M * 1 + a * n+1,n M * 2 − M * 2 a * n,n−1 − γ n (M 2 z + M 1 − M 2 a n+1,n + a n,n−1 M 2 )γ −1 n .
Observe from Theorem 3.1 that B n (z; G) is a matrix polynomial of degree at most m and A n (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.
is also a polynomial, then T = T S satisfies with G(z) = G(z) + H(z). Moreover, the matrix X n defined in (3.6), satisfies along with (3.13)-(3.16) the following relation:

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, and B n+1 (z; G) + γ −1 n B * n (z; G)γ n = (zI − α n )A * n (z; G), (3.24) where α n and β n are 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 X n in the differential equation (3.13): Under assumptions (3.2) and (3.4), with G(z) a matrix polynomial, the monic MOPRL ( P n ) n satisfy the following difference-differential relations (lowering and raising operators, respectively): P ′ n (z) + P n (z)G(z) = −B n (z; G) P n (z) + A * n (z; G)β n P n−1 (z) (3.25) and (3.26) where A n and B n are given by (3.15) and ( 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 P n , a phenomenon that has not been reported before in the theory of MOPRL. where H is defined in (3.20).
In particular, provided A n (z; G) is invertible for all z ∈ C, P n (z)H(z) + B n (z; H) P n (z) − A * n (z; H)A − * n (z; G) × P ′ n (z) + P n (z)G(z) + B n (z; G) P n (z) = 0. 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: In the scalar case (N = 1) the only smooth and unitary function s(x) on R has the form s(z) = e ip(x) , with p real-valued, and s ′ (z)s(z) = ip(z), so that the assumption that p is 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 ( P n ) n . As a consequence of the Proposition 3.6 we have the following  2) and (3.4), with G(z) a polynomial, the MOPRL ( P n ) n satisfy the following second-order differential equation provided that the inverse of A n (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 M n P ′ n , N n P n and M n P n G. In some cases, under additional assumptions on the weight, (3.30) can be reduced further, as we will see in the following section.

Illustrative examples
In this section we study a number of examples of weights W which 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 where q is a scalar real-valued function, so that with the notation (3.2) and (3.4), we may take We are interested in the case when G is a matrix polynomial. Hence, keeping up with previous hypotheses, we will assume that q is a (scalar) polynomial of even degree with real coefficients and a positive leading coefficient, and that U ′ (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 q to either 2 (the Hermite case) or 4 (the Freud case), and the degree of U ′ (x)U −1 (x) to at most 1. We start by considering the case when Taking into account the linearity (3.22) of the coefficients A n and B n in (3.15), (3.16), we can obtain the differential equation (3.13) for the general case of U ′ (x)U −1 (x) = A + 2Bx. However, finding the corresponding orthogonality weight is more involved: when A and B do not commute, solving U ′ (x)U −1 (x) = A + 2Bx is 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.

The Hermite case
For q(x) = x 2 /2, let us consider two cases, first U ′ (x)U −1 (x) = A and then U ′ (x)U −1 (x) = 2Bx. We end up this Section by discussing briefly the case of U ′ (x)U −1 (x) = A + 2Bx.
The differential equation U ′ (x)U −1 (x) = A, so that U (x) = e Ax , T (x) = e −x 2 /2 e Ax , and the weight matrix (4.1) is given by The matrix G has the form G(x) = −xI + A, and according to (3.19), The compatibility conditions of Proposition 3.5 yield and α n = 1 2 (A + γ −1 n A * γ n ). (4.5) The lowering and raising operators from Proposition 3.6 are reduced now to P ′ n (x) + P n (x)A − A P n (x) = 2β n P n−1 (x), (4.6) and Summing up the telescopic relation in (4.4) (or comparing the O(x n−1 ) term in (4.6)) gives β n = 1 2 (nI + a n,n−1 A − Aa n,n−1 ). (4.7) With the notation of Proposition 3.10, so that the differential equation (3.30) for the monic polynomials P n , orthogonal with respect to the weight (4.3), boils down to 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 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 χ = iaI for 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 where E ij is a matrix with 1 at entry (i, j) and 0 elsewhere, and a diagonal matrix For the first non-trivial example we assume that A = L and χ = iJ, so that ad A (χ) = −A (4.12) and ad 2 A (χ) = 0. Thus, H(x) = e Ax χe −Ax = i(J − Ax), and by (3.19), B n (x; H) = i(Ax − J + a n,n−1 A − Aa n,n−1 ) = i(Ax − J + 2β n − nI).

From the compatibility conditions of Proposition 3.5 we get
Jα n − α n J + α n = A + 1 2 A 2 α n − α n A 2 and where we have used relations (4.4), (4.5), (4.7) and 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 respectively, which yields by (3.29) the first order relation (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, 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án and Grünbaum in [23].
In the second case we take where L is given in (4.10) and J in (4.11). Consequently, and In this case, As to new results, we get the lowering operator, the raising operator and the first-order differential equation (3.29): Multiplying this equation by 2 and plugging it into the second-order differential equation (4.8) gives This differential equation is not of Sturm-Liouville type considered by Durán and Grünbaum 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.
The weight matrix (4.1) is now given by where we assume that B is chosen such that all the moments exist. By (4.2), T (x) = e −x 2 /2 e Bx 2 and G(x) = (2B − I)x.
The weight matrix is an even matrix function, so that all moments of odd order vanish. As a consequence, a n,n−1 = 0 and α n = 0. Then and there will be only one compatibility condition: From (4.16) we get the explicit expression for β n via 2 I − B − γ −1 n B * γ n β n = nI + 2(a n,n−2 B − Ba n,n−2 ).
The lowering and raising operators are and respectively. This yields the following second-order differential equation 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án and Grünbaum in [23]. The computations are similar and will be omitted for the sake of brevity.

Other cases
Assume now that we have a linear combination of the previous two cases, i.e.
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 A n (x; G) and B n (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 A and B commute). In general, this solution can be given in terms of the time ordered exponential U (x) = : e x 0 (A+Bs)ds :, 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 A and B are constructed (as a linear combination of the matrices L and J in (4.10) and (4.11), respectively) such that the corresponding U (x) = e (L−v 0 J )x e v 0 J x , with v 0 any real number. [22] contains also another example of a pair of non-commuting matrices A and B such that U (x) is a matrix polynomial of degree N −1, whose coefficients can be obtained recursively in terms of A and B (see Theorem A.1 therein).
An alternative way of generating examples is by considering weight matrices W of the form (4.1) with U (x) ∈ P m such that U (0) = I and det U (x) = 1 for x ∈ R. If we denote then the coefficients A k and B k are connected by simple algebraic relations. For instance, In consequence, Let us discuss some cases, depending on the degree m of the matrix polynomial U (x). In the first non-trivial example, when m = 1, we have U ( The condition U (x)U −1 (x) = I yields the following algebraic relations between the coefficients A 1 and A 2 : Algebraic manipulations of these two equations give additionally A 2 A 1 A 2 = 0 and A 2 A 2 1 = A 2 1 A 2 . Therefore The simplest case when G(x) is a matrix polynomial of degree one gives another algebraic equation, A 2 A 1 = 0. This immediately implies that A 2 2 = 0 and A 4 1 = 0, i.e. once again resulting in very strong algebraic conditions on the coefficients. In order to find an interesting example (when ad A 1 (A 2 ) = 0) we have to go at least to the dimension N = 4; for instance, Considering higher degree examples will give more and more algebraic relations among the coefficients A k , which in general are difficult to solve.
In particular, in the examples that have appeared in the literature so far, U (x) is in general a matrix polynomial of degree depending on the size N .

The Freud case
Now we have q(x) = x 4 2 , so the weight matrix is given by Again, we will consider two cases, first U ′ (x)U −1 (x) = A and then U ′ (x)U −1 (x) = 2Bx.
An expression of the second-order differential equation can be given using Proposition 3.1, but in this case it is necessary to compute the inverse of A n (x; G), which is a matrix polynomial of degree 2 with a rational function determinant in general.
Again, there is only one compatibility condition (3.23): with the lowering and raising operators P ′ n (x) + 2x P n (x)B − B P n (x) + 4xβ n P n (x) = 4 x 2 I + β n+1 + β n − 2 B + γ −1 n B * γ n β n P n−1 (x), and P ′ n (x) + 2x P n (x)B − B P n (x) = 4x 3 I + 2 2β n+1 − B − γ −1 n B * γ n x P n (x) + −4 x 2 I + β n+1 + β n + 2 B + γ −1 n B * γ n P n+1 (x), respectively. Finally, considering the O(x n−2 ) term in the lowering operator we obtain nI + 2(a n,n−2 B − Ba n,n−2 ) = 4 β n β n−1 + β 2 n + β n+1 β n − 2 B + γ −1 n B * γ n β n . This paper extends to the matrix case the methodology of derivation of differential relations for MOPRL using a RH approach. We obtain explicit formulas for the ladder operators of very general weight matrices factorized in the form W = T T * . We have also shown that in the matrix case there is an extra freedom absent in the scalar situation, which allows us to obtain a family of ladder operators, some of them of the 0-th order, which does not occur in the scalar case. Furthermore, by combining appropriately the family of ladder operators we get a family of second order differential equations satisfied by the MOPRL, some of them of the 0-th or first order (under some invertibility assumptions). This yields some new results even in the particular cases studied before in the literature. In order to keep the size of this paper reasonable, we have restricted our attention to the weights supported on the whole real line. When supp(W ) = [0, +∞) or supp(W ) = [−1, 1], we can follow the ideas exposed in Section 3.1, except that now we have to assume that either zG(z) (for the case of the semi-axis) or (1 − z 2 )G(z) (for the case of the finite interval) are matrix polynomials. For the discussion of these situations in the scalar case, where the differential equations for the Laguerre and Jacobi polynomials are obtained, the interested reader is referred to [40,Chapter 22].