A new commutativity property of exceptional orthogonal polynomials

We exhibit three examples showing that the"time-and-band limiting"commutative property found and exploited by D. Slepian, H. Landau and H. Pollak at Bell Labs in the 1960's, and independently by M. Mehta and later by C. Tracy and H. Widom in Random matrix theory, holds for exceptional orthogonal polynomials. The property in question is the existence of local operators with simple spectrum that commute with naturally appearing global ones. We illustrate numerically the advantage of having such a local operator.


A brief historical introduction
The interest of one of us in what is now called the "bispectral problem" posed and solved in [16] arose from an effort to understand and extend a mathematical miracle, uncovered by D. Slepian, H. Landau and H. Pollak at Bell Labs back in the 1960's, [43,44,55,56,57,58,59]. It is of obvious importance in signal processing and it was motivated by work of C. Shannon, see [53]. This is mentioned in the introduction to [16] and is recalled in the next few lines.
For an unknown signal f (t) supported in [−T, T ] one observes its Fourier transform F f (k) for frequencies k in the band [−W, W]. The numerically stable reconstruction of f from this data leads to the study of an integral operator in L 2 (−T, T ) with kernel given by K(t, s) = sin W(t − s)/(t − s). One needs to compute numerically many of its eigenfunctions, and this gives a very ill-conditioned problem because the eigenvalues are (except for about 4T W of them) very close together. The miracle in question is that one can exhibit a selfadjoint second order differential operator that commutes with the integral operator, has a simple and well spread out spectrum resulting in a common set of eigenfunctions. One has replaced a very ill-posed problem by a very well posed one.
The potential applications of these polynomials are quite varied: they provide a vast extension of the classical orthogonal polynomials of Jacobi, Laguerre and Hermite which feature in countless areas of mathematics, both pure and applied. They each give a basis of polynomials that are joint eigenfunctions of a fixed differential operator L of order two Lp n (x) = λ n p n (x) The main difference with the classical ones is that the index n, that indicates the degree of p n , needs not run over the entire set 0, 1, 2, 3, .... A way to increase the chances that this larger class of polynomials will be widely used is to study the extend to which they share certain properties with their classical counterparts. Such an effort will necessarily develop in an exploratory fashion. This is the spirit of this paper.
The hope at the time of [16] was that situations exhibiting the highly unusual "bispectral property", to be defined below, would give rise to extensions of the commutativity miracle mentioned above beyond the Fourier case exploited by the Bell Labs group. One says that one has a "bispectral situation" when a differential operator L has eigenfunctions f (x, k) that satisfy a differential equation in the spectral parameter k, i.e. we have (1) Lf (x, k) = Λ(k)f (x, k), as well as where B is a differential or difference operator acting on k. The spectral parameter k runs over a discrete or continuous set.
Several tools were used in [16] to classify all the situations when L has order two and B is a differential operator of arbitrary order. They include, among other things, the so called Darboux processes, the so called ad-conditions, as well as a careful study of the monodromy properties of L. The main surprise was the observation that the rational solutions of the Korteweg-deVries equation play a important role in half of the cases found in [16]. The role of its master symmetries was observed in a later paper [65].
The case when k is a discrete variable and B is a second order difference operator, as in the case of the classical orthogonal polynomials, was considered in a series of papers [27,28,29,31], where the Darboux process was applied to either the semi-infinite or doubly infinite banded matrix B to obtain the Krall polynomials or functions that extend them. Here the role of the KdV equation is taken up by nonlinear evolutions such as the Toda flows. The appearance of these integrable isospectral systems when B is a differential or a difference operator came about in these papers by using the full power of the Darboux process (see also [38,40]). Two earlier papers by M. Reach [50,51], based on his UC Berkeley thesis 1987, deal with applying the Darboux process to the second order operator L, as in [16], but allowing B to be a difference recursion of arbitrary order. Just as in [16] this leads to increasing orders in the recursion relation given by B. See also [22,41] where one has a similar situation. This is not the place to review in detail the developments just mentioned, and we just recall the very important contributions by G. Wilson, [63,64], A. Kasman and M. Rothstein, [42], as well as those of B. Bakalov, E. Horozov and M. Yakimov, [2,3,4]. More references can be found in [25,30,39].
It is clear that a certain set of common tools were used both in the study of the bispectral problem as well as in the study of exceptional orthogonal polynomials. We just mention a few of these. Apparently the first explicit mention of the bispectral property in the second area is made in [52] (see section 4 as well as the acknowledgments in this reference). It appears that the first paper to exploit the Darboux process in connection with exceptional orthogonal polynomials is [49]. The role of trivial monodromy has appeared in later papers such as [19,20]. One can see points of contact between the considerations in [22,41] and those in the bispectral problem as mentioned by these authors (see also [20]). There are many differences between these two topics: all considerations in the bispectral problem are of a local nature, whereas in the case of exceptional orthogonal polynomials the issue of the completeness of these polynomials in some appropriate Hilbert space is of importance.
Making heavy use of the bispectral property, the paper [24] establishes that the commutativity phenomenon alluded to above holds for the classical orthogonal polynomials. The paper [26] shows that the same property holds in connection with the "even family" in [16], the one connected with the master symmetries of the KdV equation. A general strategy to connect bispectrality and the commutativity property is given in [32]. A much more recent push in this direction is in [7,8,9,10].
The trivial cases of bispectrality when both physical and frequency space are the real line are given by Bessel and Airy operators (the Bessel case includes Fourier analysis). They have featured in important problems of mathematical physics, such as potential theory, electromagnetism and optics for a very long time. Both cases lead to integral operators admitting a commuting differential operator. For the Bessel case this was proved by D. Slepian, see [57]. The Airy case was observed by C. Tracy and H. Widom in the context of Random matrix theory, see [61]. In [9] one considers deformations of the Airy integral operator which preserve the commutativity property. For the Bessel case see [26]. For very recent numerical work on the eigenfunctions of the Airy integral operator which exploits the existence of the commuting one see [54]. This work has applications both in Random matrix theory as well as in optics.
In view of these more recent papers, and once one notices that exceptional orthogonal polynomials give a bispectral situation and are connected to the classical ones by applications of the Darboux process the results in the present paper are not totally unexpected. However, the nature of the explicit results given here makes it worthwhile presenting them separately. One should mention that there are matrix valued versions of the commutativity property too, see [11,12,33,34,35].
A last historical remark: many ideas appear again and again in the development of any area of mathematics and the issue of giving appropriate credit takes complicated turns. It is nice to be able to single out a paper by V. Bargmann, [1], where several of the concerns and issues in modern spectral theory originate. The very nice paper by C. Quesne [49] starts by referring to that pioneering paper. In fact, the introduction to her paper offers a rather instructive view of a rich area where many tools appear over and over again. A very complete discussion of the work of V.Bargmann, M.G. Krein, V.A. Marchenko, I.M. Gelfand and B.M. Levitan (among others) is given in the book [13].
In conclusion, we mention that the phenomenon mentioned above (in its original Fourier version) has found a rather unexpected use in a series of papers by A. Connes and collaborators in connection with the Riemann zeta function. For the most recent push in this direction, see [14]. For a commentary on this paper see [37]. For a recent use of this commutativity property (again in its original Fourier form) in connection with the Bethe ansatz and entanglement, see [5,6,15] and its references. It is interesting that many of the kernels that appear in connection with the bispectral problem play an important role in Random Matrix Theory, see [61,62]. For deformations of these kernels that preserve both phenomena see [7,8,9,10]. These observations argue for approaching the relatively new area of exceptional orthogonal polynomials by casting a wide net.

The contents of the paper
The commutativity phenomenon alluded to above was extended to the case of classical orthogonal polynomials defined on a finite set, see [47,48]. It was later observed by R. Perline that the explicit form of the commuting operator could be given a simple and unified form, see [46].
It was subsequently seen, see [33,36] that this simple form of Perline applies unchanged to other situations. In [36] one observes that (as mentioned by R. Perline) the simple form of the commuting operator may result in one that does not have simple spectrum. This happens for the so called Bannai-Ito polynomials. In [36] one shows that a more complicated form of the commuting operator, still built with an appropriate extension of Perline's construction, yields one with simple spectrum.
Since the examples of exceptional orthogonal polynomials that we consider here, namely Jacobi, Laguerre and Hermite, involve recursion relations of orders five and seven respectively, the simple form put forward by R. Perline has to be modified. This point will be illustrated in some of the examples below.
We are now in a position to describe the contents of the present paper. In section 3 we recall the definitions of the operator of time-band-time limiting as well as the operator of band-time-band limiting. In section 4 we consider the case of the exceptional Hermite polynomials which do not depend on free parameters. In section 5 we consider the simplest instance of the exceptional Jacobi polynomials. These ones depend on the usual parameters α and β and some of our results are given in terms of these parameters. Some results, which have been checked for multiple values of the parameters, are illustrated by a specific (but arbitrary) choice of them. In section 6 we consider an instance of exceptional Laguerre polynomials. Finally, in section 7, in the spirit of [11, Section 6], we conclude by exploiting the numerical pay-off of the results of the previous sections.

The operators of time and band limiting
We start with a very general setup of the time-band limiting problem for orthogonal polynomials (see for instance [24]) which will be applied later on to the different examples of sequences of exceptional orthogonal polynomials discussed in this paper. x n w(x)dx, n ≥ 0, are finite. Let (p n (x)) n≥0 be a sequence of real valued orthonormal polynomials with respect to the weight w(x). Since we will be dealing with exceptional orthogonal polynomials, we do not assume that deg p n = n. Consider the following two Hilbert spaces: The space is an isometry. If the polynomials are dense in L 2 (w), this map is unitary with the inverse F −1 : We denote our map by F to remind the reader of the usual Fourier transform. Here N 0 takes up the role of "frequency space" and the interval (a, b) the role of "physical space".
The band limiting operator, at level N acts on ℓ 2 (N 0 ) by simply setting equal to zero all the components with index larger than N . We denote it by χ N . The time limiting operator , at level T , acts on L 2 (w) by multiplication by the characteristic function of the interval (a, T ], T ≤ b. This operator will be denoted by χ T .
Consider the problem of determining a function f from the following data: f has support on the finite set {0, . . . , N } and its Fourier transform Ff is known on the set (a, T ]. This can be formalized as follows

We can combine the two equations into
To analyze this problem we need to compute the singular vectors (and values) of the operator E : ℓ 2 (N 0 ) −→ L 2 (w). These are given by the eigenvectors of the operators The operator E * E, acting in ℓ 2 (N 0 ) is just a finite dimensional matrix M , and each entry is given by The second operator S 2 = EE * acts in L 2 ((a, T ), w(x)dx) by means of the integral kernel Consider now the problem of finding the eigenfunctions of E * E and/or EE * . For arbitrary N and T there is no hope of doing this analytically, and one has to resort to numerical methods. This is a remarkably ill-conditioned problem since most of the eigenvalues are crowded together. Of all the strategies one can dream of for handling this problem, none sounds so appealing as that of finding a differential operator with simple-and spread out-spectrum which would have the same eigenfunctions as the original operators. This is exactly what Mehta as well as Slepian, Landau and Pollak did when dealing with the real line and the actual Fourier transform. They discovered (the analog of) the following properties: • For each N , T there exists a symmetric matrix L with a small number of diagonals, with simple spectrum, commuting with M . • For each N , T there exists a selfadjoint differential operator D, with simple spectrum, commuting with the integral operator S 2 = EE * .
In this paper we will see instances of exceptional orthogonal polynomials where this phenomenon holds. Once more we will see that the "bispectral property", first considered in [16], guarantees the commutativity of these two operators, a global and a local one.
For an up-to-date treatment of the important issue of computing numerically the eigenfunctions of D, see [45]. For the case of the DFT, see [23].

The exceptional Hermite polynomials
We consider the family of exceptional Hermite polynomials defined byĤ 0 = 1, where H n are the classical Hermite polynomials given by the Rodrigues formula These polynomials can also be defined by means of a determinant (see for instance [20]): The exceptional Hermite polynomials satisfy the orthogonality relation: δ n,m , n = 1, 2.

Bispectral operators.
Let us now consider the orthonormal sequence of exceptional Hermite orthogonal polynomials given by The exceptional polynomials H n satisfy the differential equation and the recurrence relation, written explicitly in [17] for a different normalization: where the coefficients α n and β n and the function Θ(x) are given by: Here, we understand that H 1 (x) = H 2 (x) = 0 as well as . This allows, in the spirit of [16], to relate Θ(x) to the appropriate Sato's τ function.
If we write and B for the semi-infinite heptadiagonal symmetric matrix: (2)), where the expression of Θ(x) is given in (4).
For fixed values of the parameters N and T one looks for a "narrow banded" commuting matrix L . One finds that there exists an heptadiagonal matrix L commuting with M . If we consider for instance the normalization L N,N = 0 and L N,N −3 = 1, this matrix is unique.
We display here the symmetric time-band limiting matrix M of size 7 that appears in the Hermite case. Here and in what follows, I will denote the identity matrix of appropriate size. .
For N = 7 we have the following expression for the corresponding commuting matrix L One may check that this matrix has simple spectrum.
Let us write B N for the truncated matrix of size N × N made up to the first N rows and columns of B and Λ N the diagonal N × N matrix Λ N = diag(0, −6, −8, . . . , −2(N + 1)). We write (7) ad X(Y ) = XY − Y X for the usual commutator and (ad X) n (Y ) = ad X (ad X) n−1 (Y ) , n ≥ 2.
Remark 4.1. The following relation holds true for every value of N , N ≥ 1: where 0 is the zero matrix of size N × N .
A non obvious consequence of this is that if one tries to write the commuting matrix L in terms of B N and Λ N for arbitrary values of N and T , in the spirit of Perline [46], one would need to use monomials of degree higher than five. We will face a much better situation in the next section. in the spirit of M. Reach in [50,51], who was the first to consider mixed bispectral situations involving differential and difference operators.
Remark 4.3. We point out that we have made no use of the Darboux process which has always played an important role in terms of the bispectral property. For the case of exceptional Hermite polynomials discussed here we notice that the differential operator L in (6) by the relation 5. The exceptional Jacobi polynomials 5.1. Bispectrality. Let α and β be the classical real Jacobi parameters, with α, β > −1 and α = β.
The existence of a five term recursion relation was established in [52]. This can be expressed by saying that the vector whereB is the pentadiagonal matrix where its n − th row is given by (...0, 0, 0,ê n ,d n ,ĉ n ,b n ,â n , 0, 0, ...).
We have by now identified the differential operator L and the difference operator B as well as the diagonal operators Λ(k) = Λ(n) and Θ(x) = (x − b) 2 (see (1) and (2)) that give us a bispectral situation for a sequence of orthonormal exceptional Jacobi polynomials.

5.2.
The commuting matrix. We find a narrow banded matrix L that commutes with the bandtime-band limiting matrix of inner products. This is illustrated in the case of size N = 7, the band limiting parameter T = 1/3 and the Jacobi parameters α = 3 and β = 4: √ 4290 The commuting matrix L of size N is unique once we ask for L N,N = 0 and L N,N −1 = 1.
As in the previous section, one writes B N for the truncated matrix of size N × N made up to the first N rows and columns of B and Λ N the diagonal N × N matrix whose entry in the position (n, n) is equal to the eigenvalue λ n = (n − 1)(n + α + β) in the differential equation (11).
Attempting to express the matrix L above as a linear combination of very simple monomials in terms of B N and Λ N , with N = 7, one obtains This expression, of course, can be written in terms of commutators and anti-commutators as in Remark 4.1.
It is worth pointing out that a similar combination, with the same monomials as above, holds true for arbitrary size N ≥ 7. This yields a non trivial extension of the results in [46]. 5.3. The commuting differential operator. For given values of N and T , consider the integral kernel In perfect agreement with the expression of L given in (15) one can build a differential operator commuting with the integral operator with the kernel given above, by taking an appropriate linear combination in terms of the operators L = L α,β and Θ(x) = (x − b) 2 (see (12) and (14)). The linear combination involves the following operators: This operator can be written in a more explicit form, namely is the orthogonality weight of the exceptional Jacobi polynomialsp (α,β) n (x) displayed above and the functions A (x), B(x), C (x) have the form Notice that the differential operator given above in (16) is selfadjoint in L 2 ((−1, T ), ρ(x)dx).
Here P 3 (x) and P 5 (x) are polynomials in x of degree three and five respectively.
In the case of P 5 (x) one gets By adding to C (x) a constant one can make, as above, one of the coefficients of P 5 (x) vanish. Here δ 5 (N ) is up to a multiplicative constant equal to (N − 1)N (N + α + β)(N + α + β + 1), the coefficients δ i , i = 0, . . . 4, are polynomials in N of degree four (not as nice as δ 5 ), λ i , i = 0, 1, 2, 4, are polynomials in N of degree two, and γ i , i = 0, 1, 2, are constants independent of N .
Remark 5.1. Using the notation in (7) one can see that the ad conditions in [50,51] expressing the bispectral property take the following form in this case where I is the identity operator and L 1 = L − (α + β + 1) 2 4 I is a proper shift of L. The operator Θ(x) = (x − b) 2 is the same as above.
This more complicated form of the ad conditions in [50,51] arises because the eigenvalue in the Jacobi case is a quadratic function of n.
Remark 5.2. For the case of exceptional Jacobi polynomials discussed here we notice that the operator L α,β in (12) is Darboux connected to the one for classical Jacobi polynomials, namely

The exceptional Laguerre polynomials
Let α > 0 and L (α) n (x), n ≥ 0, denote the classical Laguerre polynomials orthogonal with respect to the weight e −x x α in the interval (0, +∞). We consider the sequence of exceptional Laguerre polynomialŝ L The squared norms of these polynomials are given by This expression is essentially in [21, section 6.2]). We write L (α) for the sequence of orthonormal polynomials.
Turning our attention to the commuting property, we define the N × N matrix M of truncated inner products depending on a real parameter T , whose entries are given by: e −x x α (x + α) 2 dx, n, m = 1, 2, . . . , N. For each fixed value of the parameters N and T one looks for a "narrow banded" commuting matrix L .
We exhibit the commuting matrix L of size N = 7, for the special choice of the Laguerre parameter α = 7.
As before, with a proper normalization this matrix is unique.
Remark 6.2. For the case of exceptional Laguerre polynomials discussed here we notice that the differential operator L α in (17) is Darboux connected to the one for classical Laguerre polynomials, namely by the relation

The benefit of having a commuting local matrix
The previous sections have exhibited, for a collection of examples of exceptional orthogonal polynomials, a pair of matrices. The first one, a full matrix M = M T,N , is obtained by forming the inner products of the normalized OP over a restricted range in physical space, thus implementing "time limiting" with parameter T . In the resulting N × N matrix the parameter N implements "band limiting". In each case the second matrix, denoted by L = L T,N , is a narrow banded one that commutes with the first matrix, and has simple spectrum. From a numerical point of view the entire purpose of the search for this second matrix is that it reduces the problem of computing the eigenvectors of the first one, a seriously ill-conditioned one, into a very well conditioned one.
The eigenvectors of M T,N are of paramount importance since they give the singular vectors of the signal processing problem at hand, as described in Section 1.
We display below the results of some small size numerical computations that illustrate the problem of computing the eigenvectors of a full matrix, such as M T,N , some of whose eigenvalues are very close together. In each case, we give the eigenvalues of both the full matrix M T,N and those of the narrow banded matrix L T,N . It should be clear that in the case when N is large the problems indicated below get to be much worse. Our point is that they already appear for small values of N .
We use the QR algorithm as implemented in LAPACK. In each of the three situations, Hermite, Jacobi and Laguerre, we will denote by X M the matrix of eigenvectors of M T,N (normalized and given as columns of X M . We will denote by Y L the matrix of eigenvectors of L T,N (normalized and given as columns of Y L ).
In theory the eigenvectors of M T,N should agree (up to order and signs) with those of L T,N . If we compute the matrix of inner products given by Y T L X M we expect to have the identity matrix up to some permutation and possibly some signs due to the normalization of the eigenvectors which are the columns of X M and Y L .

Conclusions.
Observe that some of the entries of these matrices Y T L X M are indeed very close to the theoretically correct values, while others are terribly off. The reason is that a few eigenvalues of the full matrix of inner products M T,N are just too close together. This produces numerical instability in the computation of the corresponding eigenvectors. On the other hand, all the eigenvalues of the commuting matrix L T,N are nicely separated and the corresponding eigenvectors can be trusted.
In summary, a good way to obtain reliable numerical values for the eigenvectors of the global matrix M is to forget about M altogether and to compute numerically the eigenvectors of L . Not only we will then be dealing with a very sparse matrix for which the QR algorithm works very fast (most of the work is avoided) but the problem is numerically very well conditioned.
The main point of this paper is to show that the miracle of Slepian, Landau and Pollak, which we described in section 1, applies to the examples of exceptional orthogonal polynomials discussed here and has important consequences.

Declarations
Conflict of interest The authors declared that they have no conflict of interest.