Matrix Valued Classical Pairs Related to Compact Gelfand Pairs of Rank One

We present a method to obtain infinitely many examples of pairs $(W,D)$ consisting of a matrix weight $W$ in one variable and a symmetric second-order differential operator $D$. The method is based on a uniform construction of matrix valued polynomials starting from compact Gelfand pairs $(G,K)$ of rank one and a suitable irreducible $K$-representation. The heart of the construction is the existence of a suitable base change $\Psi_{0}$. We analyze the base change and derive several properties. The most important one is that $\Psi_{0}$ satisfies a first-order differential equation which enables us to compute the radial part of the Casimir operator of the group $G$ as soon as we have an explicit expression for $\Psi_{0}$. The weight $W$ is also determined by $\Psi_{0}$. We provide an algorithm to calculate $\Psi_{0}$ explicitly. For the pair $(\mathrm{USp}(2n),\mathrm{USp}(2n-2)\times\mathrm{USp}(2))$ we have implemented the algorithm in GAP so that individual pairs $(W,D)$ can be calculated explicitly. Finally we classify the Gelfand pairs $(G,K)$ and the $K$-representations that yield pairs $(W,D)$ of size $2\times2$ and we provide explicit expressions for most of these cases.


Introduction
Matrix valued orthogonal polynomials (MVOPs) in one variable are generalizations of scalar valued orthogonal polynomials and they already show up in the 1940s [24,25]. Since then, MVOPs have been studied in their own right and they have been applied and studied in different fields such as scattering theory, spectral analysis and representation theory [2,6,13,14,15]. In this paper we are concerned with obtaining families of MVOPs whose members are simultaneous eigenfunctions of a symmetric second order differential operator. A family of MVOPs (with respect to this pairing) is a family (P n : n ∈ N) with P n ∈ M[x] satisfying (1) deg(P n ) = n and (2) the leading coefficient of P n is invertible, for all n ∈ N. Existence of such a family is guaranteed by application of the Gram-Schmidt process on (1, x, x 2 , . . .). Moreover, up to right multiplication by GL(C N ), a family of MVOPs is uniquely determined by the weight W .
In [8] the question was raised if there exists a matrix weight together with a differential operator of degree two that has the corresponding family of orthogonal polynomials as a family of simultaneous eigenfunctions. If N = 1 then the answer is well known, we get the classical orthogonal polynomials [3]. In fact, the algebra of differential operators that have the classical polynomials as simultaneous eigenfunctions is a polynomial algebra generated by a second order differential operator, see e.g. [28].
In general, the algebra of differential operators that act on the M-valued polynomials is End(M)[x, ∂ x ]. We identify End(M) = M ⊗ M such that a simple tensor A ⊗ B acts on an element C ∈ M via (A ⊗ B)C = ABC * . A polynomial P ∈ M[x] is an eigenfunction of a differential operator D ∈ End(M)[x, ∂ x ] if there exists an element Λ ∈ M such that DP = P Λ. This is justified by the fact that we consider M[x] as right pre-Hilbert M-module.
A pair (W, D) consisting of a matrix weight and an element D ∈ End(M)[x, ∂ x ] of order two that is symmetric with respect to ·, · W is called a matrix valued classical pair (MVCP). Any family of MVOPs is automatically a family of simultaneous eigenfunctions.
Consider  [8]. After publication of [8] examples of MVCPs (W, D) with D ∈ (M ⊗ C)[x, ∂ x ] came about, arising from analysis on compact homogeneous spaces in a series of papers starting in [16] and ending in [29].
A uniform construction of MVCPs arising from the representation theory of compact Lie groups is presented in [19,32] and was inspired by [21,22,23,39]. The input datum is a compact Gelfand pair (G, K) of rank one and a certain face F of the Weyl chamber of K. For each µ ∈ F , the output is an orthogonal family of M µ = End(C Nµ ) valued functions (Ψ µ d : d ∈ N) on the circle S 1 and a commutative algebra of differential operators D(µ) for which the functions Ψ µ d are simultaneous eigenfunctions. Here, N µ is a natural number that depends on the weight µ ∈ F . Moreover the functions Ψ µ d are determined by this property and a normalization. We call Ψ µ d the full spherical function of type µ and degree d. The datum (G, K, F ) for which this construction applies is called a multiplicity free system and they are classified in [19].
We obtain families of MVOPs by multiplying the functions Ψ µ d from the left with the inverse of Ψ µ 0 . The matrix weight W µ is expressed in terms of Ψ µ 0 , some data from the irreducible K-representation and a scalar Jacobi weight that is associated to the symmetric space G/K. Conjugating the elements of D(µ) with Ψ µ 0 yields a commutative algebra D µ of differential operators for which the MVOPs are simultaneous eigenfunctions. In fact, the MVOPs are determined by this property and a normalization. The commutative algebra D µ is contained in ( The exact relation between W µ and D µ is not yet understood on the level of the polynomials, i.e. it is not clear what exactly characterizes D µ . From this point of view it is not clear what the right notion of a matrix valued classical pair should be. In the spirit of classical orthogonal polynomials the MVOPs should be determined as eigenfunctions of a commutative algebra of differential operators. The algebras D(µ) and D µ are isomorphic by definition and the first algebra is studied in [7,Ch. 9] and [27]. It would be interesting to determine its generators in our situation, where the branching is multiplicity free and the rank is one.
Since we do not have a precise description of the algebra D(µ), we content to stick to the original definition of a MVCP and we provide a method to find infinitely many examples of them. To this end we exploit the existence of a special differential operator Ω ∈ D(µ), the (image of the) second order Casimir operator Ω on G. After conjugation with Ψ µ 0 we find a second order differential operator D µ ∈ D µ that is symmetric with respect to the pairing and thus has the MVOPs as simultaneous eigenfunctions.
We show that our method is effective by applying it to low dimensional examples. First we classify the data for which we obtain MVCPs of size 2 × 2 from [19]. The short list that we obtain contains old and new items. Among the new ones are the symplectic (symmetric) pairs (USp(2n), USp(2n − 2) × USp(2)) and the spherical (but non-symmetric) pair (G 2 , SU (3)). For all these examples we provide the corresponding MVCPs below. This paper is organized as follows: In Section 2 we review the construction of families of MVOPs based on the representation theory of compact Gelfand pairs (G, K). We restrict the functions Ψ µ d to the Cartan circle A ⊂ G and we identify A = S 1 . The coordinate on S 1 is the fundamental zonal spherical function φ, normalized by two constants c, d so that , which is a function on [0, 1]. The MVOPs are basically a reflection on the three term recurrence relations of the functions Ψ µ d and Ψ µ d , introduced in this section. More precisely let Q d be defined by . Then it follows from the three term recurrence relation for Ψ µ d that Q µ d is a polynomial of degree d with non singular leading coefficient.
Section 3 is the heart of this paper. Here we discuss how the algebra of differential operators D µ comes about. By means of the bispectral property, we prove that there exist constant matrices R and S such that Ψ µ 0 satisfies the first order differential equation Remark 1.1. If we let ∂ x act on both sides (1), we obtain an instance of Tirao's matrix valued differential equation [35]. However, the techniques in [35] do not apply directly because S might have eigenvalues in −N 0 . This is indeed the case for all the examples considered in this paper.
Remark 1.2. Observe that (1), can be seen as a differential operator acting on the right on Ψ µ 0 . On the other hand Ψ µ 0 is also an eigenfunction of the radial part of the Casimir operator Ω of G acting on the left.
In Corollary 3.6, we exploit again the bispectral property of the functions Ψ µ d to deduce that the image (radial part) D µ ∈ D µ of the Casimir operator of G can be expressed in terms of R, S and an additional constant matrix coming from (G, K). More precisely we show that the polynomials where M, m are the maximum and minimum of φ| S 1 , p the period of φ and r a scaling factor. This data is provided in Table 2 for the various cases. The diagonal matrix Λ 0 is the eigenvalue of Ψ µ 0 as an eigenfunction of the Casimir operator Ω and it can be calculated for each pair (G, K). It follows that the explicit knowledge of the function Ψ µ 0 implies explicit knowledge of the corresponding pair (W µ , D µ ).
In Section 4 we discuss an algorithm to obtain an explicit expression for Ψ µ 0 . This algorithm can be implemented in in GAP [12,33] for each specific pair (G, K). Once we have a formula for Ψ µ 0 we can calculate the corresponding MVCP by differentiation and matrix multiplication. We also propose a method of finding families of MVCPs.
(1) Take a family (G n , K n , µ n ) n∈N for which the construction applies. For instance, take a constant family of Gelfand pairs and consider (G, K, nµ), where µ ∈ F , or take a canonical element of a face F , for instance the first fundamental weight, and let the Gelfand pairs vary. (2) Calculate the first so many functions Ψ µ 0 of the family until a patterns shows up. This provides an ansatz for a family of MVCPs.
(3) Show that every pair is indeed a MVCP. This is not difficult, one needs to check three equations [10,Thm. 3.1]. It turns out that in many cases the group parameter, which is a priori discrete, may vary continuously within a certain range.
In Section 5 we discuss the implementation in GAP [12] for the compact Gelfand pair (USp(2n), USp(2n − 2) × USp(2)) and the appropriate faces F . We discuss a branching rule that is necessary to implement the algorithm. The branching laws for the symplectic groups are more difficult than those for the special unitary and orthogonal groups. Indeed, the multiplicities are not only determined by interlacing condition, but also by an alternating sum of partition functions. At this point it is important to select the right irreducible K-types. Selecting an irreducible K-representation of highest weight µ ∈ F , where (G, K, F ) is a multiplicity free system, guarantees that the branching rules simplify and that the involved algebras are commutative. In fact, the whole construction of MVOPs would not work for more general irreducible K-representations.
There are two families of MVCPs of size 2 × 2 related to (USp(2n), USp(2n− 2)× USp(2)). We calculate the first family by hand. For the other family we calculated the corresponding family of MVCPs using our method. We apply the machinery once more to give an example of size 3 × 3 that is associated to this Gelfand pair.
In Section 6 we classify all possible triples (G, K, µ) that give rise to 2×2 MVCPs according to the uniform construction described in [19]. Subsequently we determine the corresponding functions Ψ µ 0 .
To indicate that our method is effective we display the MVCPs of size 2 × 2 below. Some of these MVCPs were already known (a1,b,d, see [29,31,36]) but they were obtained by different means. The other MVOPs (a2,c1,c2,g) are new as far as we know. The matrix weights are of the form We provide the expressions for Ψ µ 0 (x) and T µ instead of working out this multiplication, because the expressions become rather lengthy. The parameters n, i, m are a priori all integers for which we give the bounds in each case, see Remark 1.3. Case a: The pair (G, K) = (SU(n + 1), U (n)). We have α = n − 1, β = 0, n ≥ 2, 1 ≤ i ≤ n − 1 and m ∈ Z. We have two families of MVCPs associated to this example, depending on the sign of m, but in either case Case a1: For m ≥ 0 : .
For case c2, we do not provide a formal proof that the family of MVCPs that are produced in this way, are the ones associated to the Lie theoretical datum. For this we have to study the various representations, which is quite involved.
Case g: The pair (G, K) = (G 2 , SU(3)). There is a single 2 × 2 MVCP associated to this pair. We have α = β = 2 and We omit case d, (SO(2n), SO(2n − 1)), as it is similar to case b. We make a few remarks concerning these examples.  (15). These expressions are meromorphic in the parameters, so they remain valid. Remark 1.4. In each case the determinant of the weight is a product of powers of x and (1 − x) times a constant. On the one hand this is quite remarkable, for the weight matrices do not seem to have much structure. However, it turns out that all the weights that we construct have this property. This follows from Corollary 3.4, which also settles and earlier conjecture [32, 1.5.3].
Remark 1.5. The matrix weights W may have symmetries, i.e. they may be conjugated by a constant matrix into a diagonal matrix weight. We check whether this occurs by looking at the commutant of W µ . It turns out that we only have non trivial commutant in cases b1 and b2 for specific parameters.
Acknowledgement. The research for this paper was partly conducted when the first author visited the University of Córdoba in August and September 2012. We would like to thank the Mathematics Departments of the Universities of Nijmegen and Córdoba for their generous supports that made this visit possible.

Lie theoretical background
Let (G, K) be a pair of compact connected Lie groups from Table 1 and let g, k denote their Lie algebras. The quotient G/K is a two-point-homogeneous space (c.f. [40]) which implies that K acts transitively on the unit sphere in T eK G/K. We denote T eK G/K = p and fix a one dimensional abelian subspace a ⊂ p. The one dimensional subspace a ⊂ g is the Lie algebra of a subtorus A ⊂ G and it follows that we have a decomposition G = KAK.
Let M = Z K (A) denote the centralizer of A in K with Lie algebra m ⊂ k. Let T K ⊂ K be a maximal torus with Lie algebra t K . Then M ∩ T K = T M and AT M are maximal tori of M and G respectively, the Lie algebra of T M is denoted by t M . The complexifications of the Lie algebras are denoted by g c , . . . The (restricted) roots of the pairs (g c , a c ⊕ t M ), (g c , a c ), (m c , t M,c ) and (k c , t K,c ) are denoted by R G , R (G,A) , R M and R K respectively. We fix systems of positive roots R + G , R + (G,A) , R + M and R + K such that the natural projections R G → R (G,A) and R G → R M respect positivity.
The lattices of integral weights of G, K and M are denoted by P G , P K and P M , the cones of positive integral weights by P + G , P + K and P + M . The theorem of the highest weight implies that the equivalence classes of the irreducible representations are parametrized by the cones of positive integral weights. Given λ ∈ P + G we denote by π G λ : G → GL(V λ ) an irreducible representation of highest weight λ. The restriction π G λ | K decomposes into a finite sum of irreducible K-representations and for µ ∈ P + K we denote the multiplicity by m G,K Definition 2.1.
K be a face, i.e. the N-span of some fundamental weights of K. A triple (G, K, µ) is called a multiplicity free triple if (G, K, µ) is a multiplicity free triple for all µ ∈ F .
The notion of a multiplicity free system can be considered for any compact Lie group G with closed subgroup K. In [19,32] it is shown that for (G, K, F ) to be a multiplicity free triple, the pair (G, K) is necessarily a Gelfand pair. Furthermore, the multiplicity free systems (G, K, F ) with (G, K) a Gelfand pair of rank one are classified by the rows of Table 1. The multiplicity free systems with (G, K) a compact symmetric pair have been classified in [18]. Table 1. Compact multiplicity free systems of rank one. In the third column we have given the highest weight λ sph ∈ P + G of the fundamental zonal spherical representation in the notation for root systems of Knapp [20, App. C], except for the case Let (G, K, F ) be a multiplicity free system from Table 1, let µ ∈ F and define P + The structure of the set P + G (λ) is important for the construction of the familie of matrix valued orthogonal polynomials that we associate to (G, K, µ). In the special case µ = 0 we have P + G (0) = Nλ sph , where λ sph is called the fundamental spherical weight. We have indicated the weights λ sph in Table 1.
The complexified Lie algebra g c has a decomposition g c = k c ⊕ a c ⊕ n + , where n + is the sum of the root spaces of the positive restricted roots R + (G,A) . If (G, K) is symmetric this is just the Iwasawa decomposition. For the two non-symmetric cases see e.g. [32].
, see e.g. [19]. It follows that the natural projection q : P G → P M induces a map P + G (µ) → P + M (µ) which turns out to be surjective, see [19]. On the other hand, if λ ∈ P + G (µ) then λ + λ sph ∈ P + G (µ), which follows from an application of the Borel-Weil Theorem. Define the degree d : is called the µ-well and B(µ) the bottom of the µ-well.
There is an isomorphism The proof of Theorem 2.2 is based on a case by case inspection, see [19,32]. Let µ be the partial ordering on P + G (µ) induced from the partial ordering (≤, ) on Nλ sph + B(µ).
Fix a multiplicity free system (G, K, F ) from Table 1 and fix µ ∈ F . Let π K µ : K → GL(V µ ) be an irreducible representation of highest weight µ. Let R(G) denote the (convolution) algebra of matrix coefficients of G and define the K × K-action Furthermore, note that E 0 is a polynomial algebra and that E µ is a free, finitely generated E 0 -module. In fact, E µ ∼ = E 0 ⊗ C |B(µ)| as E 0 -modules. Let λ ∈ P + G (µ) and let π G λ : G → GL(V λ ) denote the corresponding representation space. Then V λ = V µ ⊕ V ⊥ µ and we denote by b : V µ → V λ a unitary K-equivariant embedding and by b * : V λ → V µ its Hermitian adjoint.
It is clear that the elementary spherical functions have the desired transformation behaviour. We equip the space E µ with a sesqui-linear form that is linear in the second variable, with dg the normalized Haar measure. As a consequence of Schur orthogonality and the Peter-Weyl Theorem we have the following result.
Theorem 2.7. The matrix valued polynomial Q d is of degree d and has invertible leading coefficient.
Corollary 2.3 implies that the matrices A d are upper triangular and the non vanishing of c µ λ,λ+λ sph that the diagonals are non-zero.
The pairing Q, Q ′ = G Q(φ(g)) * W µ (φ)Q ′ (φ(g))dg is a matrix valued inner product, see e.g. [19]. Note that all the functions in the integrand are polynomials in φ, and φ is K-biinvariant. In view of G = KAK and the integral formulas for this decomposition, we have where x = cφ + d (with constants c, d depending on the pair (G, K)) is a normalization such that x attains all values in [0, 1]. The factor (1 − x) α x β is the ordinary Jacobi weight that is associated to the Riemann symmetric space Now we come to a different discription of the family (Q d : d ∈ N), one that allows us to transfer differentiability properties of the elementary spherical functions to similar properties of the matrix valued polynomials.
The spherical functions Φ ∈ E µ are determined by their restriction to A and consists of diagonal matrices. More precisely, with respect to a basis of the Msubrepresentations of V µ , the matrix Φ(a) is block diagonal, and every block is a constant times the identity matrix of size dim ν. Sending such a matrix to a vector containing these constant provides an isomorphism u : is the matrix valued function whose columns are the vector valued functions Ψ µ λ(d,ν) , with ν ∈ P + M (µ), and where λ : , with T µ the diagonal matrix whose entries are dim ν. Moreover, we know precisely which matrix coefficients occur in the entries of the functions Ψ µ d .
Theorem 2.9. The entries of Ψ µ d are indexed by the set P + M (µ). Hence, the entry The proof is immediate from the definition of Ψ µ n . Note that the construction of the functions Ψ µ d can also performed for the complexified pair (G C , K C ). We obtain a function on A C ∼ = C × that takes values in End(C |B(µ)| ) and we denote it with the same symbol, Ψ µ d : C × → End(C N ). For later reference, we state the following result concerning the entries of Ψ µ d , of which the proof is straight forward.
Proposition 2.10. The entries of the function Ψ µ λ : A C → C |B(µ)| are Laurent polynomials in C[z]. The maximal degree is less than or equal to |λ(H A )|, where H A is defined by A = a/2πiH A Z. In fact, the maximal degree occurs precisely in the entry labeled with λ| tM .
The functions Ψ µ d are analytic, as the entries are matrix coefficients. Moreover, they satisfy the three term recurrence relation where the matrices A d , B d , C d are the same as in (4).
The operator ∆ µ is a second order difference operator that has Ψ µ d as eigenfunction and with eigenvalue φ.

Differential properties
In this section we discuss the second order differential operator that we obtain from the quadratic Casimir operator Ω of the group G. Moreover, using the fact that Ψ µ 0 is an eigenfunction of Ω whose eigenvalue is a diagonal matrix and that the functions Ψ µ n satisfy a three term recurrence relation, we deduce that Ψ µ 0 satisfies a first order differential equation. From the singularities of this equation we deduce that Ψ µ 0 (a) is invertible whenever a ∈ A is a regular point for φ| A . In this section we work with spherical functions on the complexified Lie groups G C and A C .
Let U (g c ) be the universal enveloping algebra of g c and let U (g c ) K denote the algebra of differential operators that are invariant under the pull back of right The irreducible representations of D(µ) correspond to irreducible representations of g c whose restriction to k c has a subrepresentation of highest weight µ, see e.g. [7, Théorème 9.1.12].
The differential operators D ∈ D(µ) leave the space of µ-spherical functions invariant. As the µ-spherical functions are determined by their values on A, in view of G = KAK and the transformation behavior of the µ-spherical functions, Theorem 3.1. For every element D ∈ D R (µ) and every d ∈ N there is an element is a representation. The family of representations {Λ d } d∈N separates the points of the algebra D(µ).
Proof. The first part is proved in [19,32]. The matrix Λ n (D) is diagonal and its entries are Λ n (D) ν,ν = π G λ(n,ν) (D). For the second statement, suppose the converse. Then there are two differential operators D, D ′ that have the same eigenvalues and it follows that D − D ′ acts as zero on the elementary spherical functions. This implies that D − D ′ = 0.
It follows that for any D ∈ D R (µ) and ∆ µ (defined in (6)), the triple (∆ µ , D, {Ψ µ d : d ∈ N}) has a bispectral property, i.e. the operators ∆ µ and D have the members of the family {Ψ µ d : d ∈ N} as simultaneous eigenfunctions (albeit in different variables). We will use this interplay in Section 3. For more on the bispectral property, see e.g. [9,17]. In the case µ = 0 we have |B(µ)| = 1 and we denote the eigenvalues by lower case letters, Dφ = λ(D)φ.
The Casimir operator on G is given as follows. Let {X i : i = 1, . . . , dim g c } be a basis of g c and let { X i : i = 1, . . . , dim g c } be a dual basis with respect to the Killing form κ on g c . Then Ω = i,j κ(X i , X j ) X i X j . On the g c -representation V λ the Casimir operator Ω acts with the scalar κ(λ, λ) + 2κ(λ, ρ G ), where ρ G is half the sum of the positive roots of G. The image of Ω in D(µ) is denoted by Ω or by Ω(µ) if we want to indicate on which function space we let it act. Then where c(z) is a meromorphic function on A C , F µ (z) is a meromorphic function with matrix coefficients and r is a constant that depends on the pair (G, K). For the particular case µ = 0 we have F 0 (z) = 0. For the symmetric pairs these statements follow from [41, Proposition 9.1.2.11]. For the two non-symmetric pairs see [19] or [32, 3.4.28]. In what follows we only consider the operator R(µ, Ω), the radial part of the image of the second order Casimir operator in D(µ). The eigenvalues are denoted by Λ d for general µ ∈ F and λ d for µ = 0, i.e. we have Theorem 3.3. The function Ψ µ 0 satisfies the first order differential equation Proof. We need the identities (6), (8) and (9) for d = 1, 2 to prove this result. Let the operator r(z∂ z ) 2 act on both sides of the equality φ(z)Ψ µ 0 (z) = Ψ µ 1 (z)A 0 + Ψ µ 0 (z)B 0 and work out the differentiation. The derivatives of order > 1 can be written in terms of φ, Ψ µ 0 and Ψ µ 1 and their first order derivatives. We get Equating and using the three term relation and its derivative, we find . Using the three term recurrence relation once more we get . Plugging in R and S yields the desired equation.
Note that the matrix R measures the fact that the matrices A n of the recurrence relations are diagonal or not. In the symplectic case the recurrence matrices are not diagonal in general, see Sections 5.2 and 5.3. For the orthogonal groups (SO(m + 1), SO(m)) the recurrence matrices A n are diagonal. For m = 3 this follows from [22,Thm. 3.1] or [30,Thm. 9.4]. In general this follows from the description of P + G (µ) and the decomposition of the tensor product of a general irreducible representation and the fundamental spherical representation, see [32,Ch. 2.4]. Let Finally let A C,µ−reg denote the set of points where det Ψ 0 (z) = 0. Since the weight function W µ is polynomial in φ and moreover, the highest degree polynomial occurs precisely once in each column and once in each row, the determinant of W µ is a polynomial in φ of positive degree. It follows that A C,µ−reg ⊂ A C is a dense open subset. See also [32]. Proof. The linear system (10) has singularities precisely in A C \A C,reg . Hence, locally in A C,reg , there exists a fundamental solution matrix which is holomorphic and invertible. By Theorem 3.3 the function Ψ t 0 has the same properties on the set A C,µ−reg , whose intersection with A C,reg is dense in A C . The claim follows.
We can now conjugate the image Ω(µ) ∈ D(µ) of the quadratic Casimir operator Ω with the function Ψ 0 to obtain a differential operator of order two for the family of matrix valued orthogonal polynomials. Since Ω(µ) has real eigenvalues the resulting operator is symmetric with respect to the pairing (5). Indeed, the matrices Q d , Q d are diagonal matrices.
Proof. It is a straightforward computation that It follows from (7) and Lemma 3.2 that the coefficient of Q(z) in (11) is exactly Λ 0 . On the other hand, the coefficient of (z∂ z Q(z)) in (11) is obtained from Theorem 3.3. Now the theorem follows by multiplying with Ψ −1 0 from the left on both sides of (11). Then Ψ µ 0 satisfies the (right)-first order differential equation Here M, m are the maximum and minimum of φ| S 1 , see Lemma 3.2 Proof. The proof is a straightforward consequence of Lemma 3.2, Theorem 3.5 and (8). In Subsection 4.3 we provide all the data r, p. The scaling r is determined by comparing the radial part of Ω(0) to the hypergeometric differential operator for the Jacobi polynomials on G/K. Basically we only need to know (α, β).
Remark 3.7. The operator D µ in (13) is a matrix valued hypergeometric equation [35]. In the scalar case, the polynomial eigenfunctions of (13) can be written in terms of hypergeometric series. In the matrix valued setting, in order to give a simple expression of Q d as matrix valued hypergeometric series it is necessary to perform a deeper analysis. See Corollary 5.6 for the case of the symplectic group.
Remark 3.8. Theorem 3.5 allows us to calculate the conjugation of the Casimir operator with the function Ψ µ 0 for induvidual cases, without calculating the radial part of the Casimir oparator. The latter is in general very technical so Theorem 3.5 makes it much easier to generate explicit examples of differential operators. For an explicit expression of the differential operator D µ we only need to know the eigenvalue λ 1 and an explicit expression of the function Ψ µ 0 . Indeed, using Lemma 3.3 we find expressions for R and S, using computer algebra, and these, together with λ 1 gives an expression for D µ by (13).
We conclude that a numerical expression of the functions Ψ µ 0 allows us, using computer algebra, to generate examples of a matrix valued classical pairs (W µ , D µ ). Remark 3.9. The operator D µ is Theorem 3.6 is given explicitly by • For SO(2n + 1):

A method to calculate MVCPs
In this section we describe an algorithm to calculate the functions Ψ µ d from Definition 2.8. We have implemented this algorithm in the computer package GAP for the symmetric pair (USp(2n), USp(2n − 2) × USp(2)) (see Subsection 5 for this pair and [33] for the GAP-files). From the description of the algorithm it is clear how to implement it for other multiplicity free triples in Table 1. In the §4.1 we discuss the general algorithm to calculate Ψ µ n . In §4.2 discuss the implementation of the algorithm in GAP. Finally we provide the necessary data from the compact Gelfand pairs to calculate the MVCPs.
Put all the entries in the j-th column 10. Build Ψ µ d with all the columns j = 1, . . . , N µ .
Step 1. Determine P + M (µ). The set P + M (µ) parametrizes the the elementary spherical functions of a fixed degree d (Theorem 2.2) and also the entries of the elementary spherical functions restricted to A (Theorem 2.9). Thus, the number of elements N in P + M (µ) determines the size of the matrices Ψ µ d (a), a ∈ A. Steps 2 and 3. Determine B d (µ). This set determines the spectrum P + G (µ) and in turn parametrizes the elementary spherical functions.
Step 5. Determine a K-equivariant embedding γ : V µ → V λj . To this end we calculate the weight spaces V λj (µ ′ ), where µ ′ is a weight in P G that projects onto µ ∈ P + K , i.e. with q(µ ′ ) = µ. If rankG = rankK then µ ′ = µ. The highest weight vector v 0 ∈ V µ is annihilated by all the simple root vectors of K and the image γ(v 0 ) is the unique vector with this property. Hence we calculate which is a one-dimensional subspace of V λj because m G,K λj (µ) = 1. Fix a non-zero element v ′ 0 in (14). Put γ(v 0 ) = v ′ 0 and define γ : V µ → V λj by stipulating that γ is K-equivariant. Letting the root vectors in k C of the negative simple roots act on v ′ 0 gives a basis of γ(V µ ) ⊂ V λj .
Step 7. Determine a vector v λj ,µ,νi in γ(V µ ) of highest M -weight ν i . This is a similar to the first part of Step 6.
Step 8. Determine the matrix coefficient e it → a t v λj ,µ,νi , v λj ,µ,νi . By Proposition 2.10 we know that any matrix coefficient of V λj restricted to A is a Laurent polynomial whose terms are of degree k with |k| ≤ λ(H A ). On the other hand we have a t = exp(iH A t), which implies It suffices to calculate H k A v λj ,µ,νi , v λj ,µ,νi V λ for k ≤ 2|λ(H A )|. Indeed, the Laurent polynomial m λ v,w | A is determined by the the first 2|λ(H A )| terms of its power series. To get an invariant inner product on the representation spaces we have implemented the Shapovalov form, see e.g. [34].
Steps 9 and 10. The matrix coefficients we obtain in Step 8 are put in a column vector of length N µ and these N µ column vectors are in turn put as columns in an N µ × N µ -matrix.

Implementation in GAP. In
Step 1 we need to determine the P + M (µ). Descriptions of P + M (µ) are given by classical branching rules for the lines 1,2 and 3, see [32]. For line 4 see Section 5, for line 5 see [1] and for the two remaining lines see [19]. Upon choosing suitable tori, the branching rules amount to interlacing conditions of strings of integers which are easily implemented in GAP.
To implement the bottom B(µ) for Steps 2 and 3 we need a precise description. For the first three lines of Table 1 see [4,32] and for the other lines see [19]. Once B d (µ) is known for d = 0, it is known for all d. The bottom depends (piecewise) affine linearly on P + M (µ) and is thus easily implemented in GAP. Having the highest weights of the involved irreducible G, K and M -representations we can do all the Lie algebra calculations in the appropriate G-representation space V λ using GAP. This settles Steps 5 and 7.
The linear system that we need to solve in Step 8 to determine the coefficients of the Laurent polynomials m λ vν i ,vν i | A is easily implemented. The actual output of the algorithm is a number of N 2 strings of real numbers, each string representing a Laurent polynomial m λ vν i ,vν i | A . A simple loop provides the implementation of Steps 9 and 10.
Remark 4.1. The representation spaces that are used soon become very large, which makes the implementation only suitable for relatively small K-types µ. In this case, small means that the calculation is performed in a reasonable time. It depends on, among other things, the dimensions of the representations spaces V λ , λ ∈ B(µ). It would be interesting to make these things more precise, once the algorithm is improved. For instance, of the spaces V λ , λ ∈ B(µ), we only need a few vectors, whereas in our implementation we first need the whole space V λ to calculate these few vectors. Unfortunately, we do not see how we can realize these kind of improvements now.

4.3.
Obtaining the MVCPs. Once we have constructed the function Ψ µ 0 whose entries are polynomials in cos(t) we follow the next steps (a) We make the change of variables x = (cos(pt)+1)/2 so that Ψ µ which gives the matrices R, S according to Theorem 3.3. (c) The differential operator D µ is now given by Remark 3.9. The missing data m, M, p is contained in the fundamental spherical function φ (see the proof of Lemma 3.2), which, as well as λ 1 , we provide in Table 2. The eigenvalue Λ 0 has to be calculated in the individual cases using for example Weyl's dimension formula. We do this in the Sections 5 and 6 for the 2 × 2-examples. (d) To obtain an expression for the weight W µ (x) = W µ (x)w(x) we need the scalar weight w(x) = (1−x) α x β and the matrix T µ that has the dimensions of the irreducible M -representations ν ∈ P + M (µ) on the diagonal, see (5) and the discussion following Definition 2.8. The parameters (α, β) of the scalar weights are provided in Table 2. The weight is given by Hence, knowledge of (Ψ µ 0 , φ, T µ , λ 1 , Λ 0 ) allows us to find explicit expressions of the pair (W µ , D µ ). (e) A formal proof that (W µ , D µ ) is a MVCP can be done showing that D µ is symmetric with respect to W µ . This boils down to prove that the following symmetry equations, [10, Theorem 3.1], hold true where We observe that the boundary conditions in [10, Theorem 3.1] are always satisfied in our case.

The symplectic case
Let n ≥ 3 and G = USp(2n), K = USp(2n−2)×USp(2). The branching rules for G to K are due to Lepowsky [26]. Let K 1 = USp(2) × USp(2n − 4) × USp(2) ⊂ K and M = USp(2) × USp(2n − 4) ⊂ K 1 , where the embedding K 1 ⊂ K is the canonical one and where M ⊂ K 1 is given by (x, y) → (x, y, x). The branching rules for K to M are due to Baldoni-Silva [1]. If we choose the K-type in a 2dimensional face then the branching rules simplify considerably. We employ the standard choices for roots and weights of the symplectic group [20, App. C]. We have rankG = rankK = rankK 1 = n and rankM = n − 1. The weight lattices of G, K, K 1 are equal to P = Z n . Let {ǫ i : i = 1, . . . , n} denote the standard basis of P . The fundamental weights of K are ω i = i j=1 ǫ i for i = 1, . . . , n−1 and ω n = ǫ n . The fundamental weights of K 1 are ξ 1 = ǫ 1 , ξ i = i j=2 ǫ i for i = 2, . . . , n − 1 and ξ n = ǫ n . The fundamental weights of M are identified with η 1 = 1 2 (ǫ 1 + ǫ n ) and η i = i j=2 ǫ j for i = 2, . . . , n − 1. The fundamental weights generate the semi group of integral dominant weights P + K , P + K1 and P + M which in turn parametrize the equivalence classes of irreducible representations of K, K 1 and M . 5.1. Branching rule from K to M . The branching rule from K to M can be calculated in two steps: first from K to K 1 , then from K 1 to M . A crucial ingredient of the branching rule is the partition function p Ξ : Z n−1 → N, where Σ = {ǫ i ± ǫ 1 : i = 2, . . . , n − 1} and Proposition 5.1. Let µ = xη i +yµ j with i < j and x, y ∈ N. Write µ = n k=1 b k ǫ k and let ν = n k=1 c k ǫ k ∈ P + K1 . Define Proof. The statement is Lepowsky's branching rule for USp(2n − 2) to USp(2) × USp(2n − 4) with the additional information that we have control over the multiplicity. The support of the function m K,K1 µ : P + K1 → N is contained in the set that is determined by condition (1) and (4) and the additional condition that n−1 k=1 C k is even, see [26,Thm. 2]. In this case the multiplicity is given by where p Ξ is the partition function (16). See [20,Thm. 9.50] for a proof of an equivalent statement (the roots are permuted because the embedding is different). If m K,K1 µ (ν) = 1 then ν is in the support of m K,K1 µ which implies that C k ≥ 0 for k = 1, . . . , n − 2. The condition on µ implies that C k = 0 unless k = i or k = j, hence (2) is satisfied. Condition (1) is trivially satisfied. It remains to check (3), which we do below. Conversely, if (1) -(4) are satisfied, then ν is in the support of m K1,K µ , because C k = 0 unless k = i or k = j. It remains to show that (3) implies m K,K1 µ (µ) = 1, which we check by calculating (17). We distinguish two cases: i = 1 and i > 1. Suppose that i = 1. Then C k = 0 unless k = 1 or k = j. Without loss of generality we may assume that j = 2 (if j = n then we take C 2 = 0) and Ξ = {ǫ 2 ± ǫ 1 }. Then Both terms are zero or one, because |Ξ| = 2. The first term is one if and only if C 1 − C 2 ≤ c 1 ≤ C 1 + C 2 and C 1 + C 2 − c 1 even. The second term is minus one if and only if c 1 ≤ C 2 − C 1 − 2. As c 1 ≥ 0 in the first place, we see that m K,K1 µ (ν) = 1 if and only if (1) -(4) are satisfied. Suppose that i > 1. Without loss of generality we assume that i = 2 and j = 3 and Ξ = {ǫ 1 ± ǫ 2 , ǫ 1 ± ǫ 3 }. Then A similar calculation shows The quantities (18,19) can only be nonzero if C 2 + C 3 − c 1 is even. In this case (18) is the number of integral points of the intersection of the line ℓ = {(B 2 , B 3 ) : (19) is equal to the number of points of the intersection of the same rectangle with the line ℓ + (1, 0). It follows that m K,K1 The branching from K 1 to K is equivalent to the branching USp(2) × USp(2) to the diagonal USp(2). Given a K 1 type (c ′ 1 ǫ 1 + c 2 ǫ 2 + · · · + c n−1 ǫ n−1 + c ′ n ǫ n ), the M -types that occur upon restricting are of the form (c 1 ǫ 1 +c 2 ǫ 2 +· · ·+c n−1 ǫ n−1 +c 1 ǫ n ), Corollary 5.2. The matrix valued orthogonal polynomials of size 2×2 are obtained from µ = ω 1 or µ = ω n−1 .
If 1 < i < n − 1 then µ = (x + y)ω i then (c i , c i+1 ) ∈ {(0, 0), (x + y, 0), (x + y, x + y)} lead to different representations of K 1 . Case 2. We have µ = xω i + yω n with xy = 0. If 1 < i < n − 1 then we find three K 1 types that occur upon restriction. If i = 1 then we can choose c 2 = x and c 2 = 0. The first choice implies c 1 = x, the second c 1 = 0. We find two K 1 -types, one of them decomposes in at least two M -types upon restriction.
Case 3. If x, y are both non zero, then we have at least three K 1 -types. Indeed, if either x or y is zero, then we are in Case 1. Assume xy = 0. Then µ, µ + ǫ 1 − ǫ i , µ + ǫ 1 − ǫ j are three K 1 types by Proposition 5.1.
We conclude that µ = ω 1 and µ = ω n−1 are the only K-types that give matrix valued orthogonal polynomials of size 2 × 2.

5.2.
The algorithm for a 2 × 2-case. We consider (USp(2n), USp(2n − 2) × USp (2)) and the irreducible USp(2n − 2) × USp(2)-representations with highest weight µ in a two dimensional face. The complexified Lie algebra of USp(2n) is denoted by sp 2n (C). We calculate the spherical function φ and Ψ µ 0 . The realizations of the fundamental representations that are involved in this calculation are described in [11]. We employ the same notation as in [11]. In particular, the root vectors of ǫ i − ǫ j and 2ǫ k are denoted by X i,j and U k respectively.
The spherical representation λ sph = ̟ 2 = ǫ 1 + ǫ 2 is realized in the kernel of  (2) is of weight zero because the groups G and K are both of semisimple rank n. We look for a vector that is killed by sp 2n−2 (C) ⊕ sp 2 (C) and it is easily seen that e 1 ∧ e n+1 + . . . + e n−1 ∧ e 2n−1 − (n − 1)e n ∧ e 2n is killed by this subalgebra.
The torus A has as its Lie algebra a = R · (X 1,n − X n,1 ). Put Then A = {a t |t ∈ [0, 2π)} and v 0 , a t v 0 = n 2 cos 2 (t) − n and it follows that φ(a t ) = (n cos 2 (t) − 1)/(n − 1). Now we calculate the function Ψ ω1 0 by means of the algorithm. Upon identifying V ̟1 ∼ = C 2n we have v ǫi = e i and v −ǫi = e n+i .
Step 4-10: we calculate Ψ ω1 ̟3 (a t ). Let V = V ̟1 as above. We realize V ̟3 in the kernel of the map ϕ 3 : In 3 V , the weight vectors of weight There are 2n(n−1) of them. The kernel of this map is of dimension 2n 3 −2n = 8 n 3 +2(n−2)n, the multiplicities of the short roots are n − 2 and there are 2n of them. It follows that the restricted map 3 V (±ǫ i ) → V has an (n − 2)-dimensional kernel. Since In order to calculate the embedding V ω1 into V ̟3 we need to calculate α∈R + K ker(e α | ker(ϕ 3 | 3 V (±ǫ i ) ) ).
Proof. The result is proved by following the steps in Subsection 4.3.
Remark 5.5. Let us consider a differential operator of the form where C, U and V are d × d matrices and F : C → C d is a (column-)vectorvalued function which is twice differentiable. It is shown by Tirao [35] that if the eigenvalues of C are not in −N, then the matrix-valued hypergeometric function 2 H 1 defined as the power series converges for |z| < 1 in M d (C). Moreover, for F 0 ∈ C d the (column-)vector-valued function is a solution to (21) which is analytic for |z| < 1, and any analytic (on |z| < 1) solution to (21) is of this form.
In the following Corollary we write the monic MVOPs with respect to W ω1 in terms of the matrix-valued hypergeometric functions 2 H 1 .
Corollary 5.6. The unique sequence of monic MVOP {P d } d≥0 with respect to W ̟1 is given by Here e i is the standard basis vector.
Proof. Since the eigenvalues of 2 + 2 S are not in −N, we can apply Remark 5.5. By looking at the leading coefficient, it is easy to see that the eigenvalue of D ω1 for P d is given by Λ Observe that the matrix Λ M d is diagonalizes after conjugation by M . Therefore the columns of P d M −1 are a solution of an equation of the form (21). Now the proof follows by writing the appropriate solution to (21). 5.3. The second family of 2 × 2 MVOP. In this subsection we use our implementation in GAP ( [33]) to obtain a one-parameter sequence of MVOP associated with the spherical functions of type µ n = ̟ n−1 . For this purpose we compute Ψ µn 0 for small values of n and we make an ansatz for the general expression of Ψ µn 0 . Finally we prove that (W µn , D µn ) is a MVCP by solving the symmetry equations (15).
. The matrices Ψ µn 0 and T µ are the building blocks of the weight matrix W µn . In the following lemma, we compute explicitly the first order differential equation for Proof. The lemma follows by computing explicitly Then S is minus the coefficient of degree 0 in x and R is the coefficient of degree 1 in x.
We can now apply Theorem 3.3 to obtain the differential operator D µn . Using the symmetry relations we (15) we check that the thus obtained pair (W µn , D µn ) is a MVCP.
Theorem 5.9. Let W µn and D µn be defined by with R and S given in Lemma 5.8. Then (W µn , D µn ) is a MVCP.
Remark 5.10. Observe that, as in the previous example, the eigenvalues of S are 0 and 1/2. This implies that the eigenvalues of 2 + 2 S are 2 and 3 so that the polynomial eigenfunctions of the differential operator in 5.9 can be written in terms of matrix valued hypergeometric functions [35].

5.4.
A family of 3 × 3 MVOP. In this subsection we fix the multiplicity free triple (USp(6), USp(4) × USp (2)) and we study the one-parameter family of Ktypes µ j = ̟ 1 + j̟ 3 . For small values of j ∈ N we can use our implementation in GAP to compute Ψ µj 0 . We use this information to make the following ansatz of Ψ µj 0 for all j: For all j ∈ N, we have As in the previous example, we obtain the matrices R and S in the differential equation (13) for Ψ µj 0 .
Now we proceed as in Section 4.3 to construct the MVCP associated to µ j .
Theorem 5.11. Let W µn and D µn be defined by Then (W µn , D µn ) is a MVCP.
Remark 5.12. The eigenvalues of S are (i + 1)/2, i/2, (i − 1)/2 so that 2 + S has always positive eigenvalues. This ensures that the hypergeometric operator (23) fits into the theory developed in [35]. However, further analysis is required to find an expression of the orthogonal polynomial with respect to W µn as matrix hypergeometric functions.
Proof. We have to prove that the indicated irreducible K-representations decompose into two irreducible M -representations upon restriction to M . In all cases this follows from the branching rules K ↓ M described in [4] (a,b,d), Section 5 (c) and [19] (g). To exclude all the other cases for the pairs (G, K) in (bf a,b,d) we refer to the same branching rules. These are all given by interlacing conditions and the statement is easily checked. Roughly speaking, the K-weight is given by a string of ordered integers and as soon as we make more than one jump, of a jump larger than one, we will find more than two M -weights that interlace. Case (iv) is dealt with in Section 5. For (g) we have to check the branching SU(3) ↓ SU(2) which can also be done in stages, SU(3) ↓ U(2) ↓ SU(2). On the first branching case (bf a) applies. We are left to show that the other multiplicity free systems are excluded. We refer to [19,1] for the branching rules G 2 ↓ SU(3) and Spin(9) ↓ Spin(7) respectively. These are the appropriate branching rules K ↓ M for the two remaining cases (G, K) in Table 1 and one readily checks the assertion.
The maps f, f ′ are M -equivariant, i is an M -isomorphism and γ, γ ′ are the embeddings of the K-representation spaces and are thus K-equivariant. Let w ν , w ν ′ denote the highest weight vectors of the M -representations. We need to calculate the four vectors x 1,1 = f (w ν ), x 2,1 = f (w ν ′ ), x 1,2 = f ′ (w ν ) and x 2,2 = f ′ (w ν ′ ) and subsequently the inner products x i,j , a t x i,j . Together with information from Table 2 we calculate the data ( Ψ µ 0 , φ, T µ , λ 1 , Λ 0 ). Now we follow the steps in Section 4.3 to obtain explicit expressions for the corresponding pair (W µ , D µ ).
In the remainder of this section we provide all the functions Ψ µ 0 that can be obtained from Table 1 using our method (Note that case c is dealt with in Section 5). Once we have all the pairs (W, D) we make some final remarks about the parameters that are involved, about the decompositions of the weights and about the question whether we can calculate the corresponding MVOPs in terms of Tirao's matrix valued hypergeometric function [35]. Case a: (G, K) = (SU(n + 1), U(n)). In this case we need to distinguish between two cases. The first case has been already investigated in a series of papers ending with [29].
Case a2: This case was not considered in the literature before as far as we know. We take the K-types µ = ̟ i + m̟ n but now we let m ∈ Z <0 . Note that although all ingredients needed for this case are already given in [29], the weight matrix constructed there does not have finite moments for negative values of m, thus it is excluded [29, Section 6]. Our case is essentially a conjugation of the case considered in [29] and does not have any problem of integrability. The function Ψ µ 0 is given by Cases b and d. We only treat (G, K) = (SO(2n + 1), SO(2n)) and µ = ω i with i ∈ {1, . . . , n − 2}. This case is essentially the same as the one in Section 5 but simpler. We use the notation of [11,Ch. 19]. The K representation is of highest weight ω i for indicated i is realized in i C 2n and this space decomposes as i−1 C 2n−2 ⊕ i C 2n−2 as M -representation, where M ∼ = SO(2n−1). It follows that P + M (ω i ) = {η i−1 , η i }. On the other hand, B(ω i ) = {̟ i , ̟ i+1 } and λ(0, η i ) = ̟ i+1 . Everything is now explicit and the maps f, f ′ of the commutative diagram (24) are easily calculated on highest weight vectors. We have A = {a t = exp(tH A ) : t ∈ R} where H A = X 1,n − X n,1 . We find Ψ ωi 0 (a t ) = cos(t) 1 1 cos(t) .