Spectra of Jacobi operators via connection coefficient matrices

TL;DR

This paper develops a method to compute the spectrum and spectral measure of Jacobi operators with compact perturbations using connection coefficient matrices, reducing the problem to polynomial root finding and establishing computability with error guarantees.

Contribution

It introduces an explicit polynomial-based approach for spectral computation of perturbed Jacobi operators and proves the spectrum's computability within the Solvability Complexity Index framework.

Findings

Spectrum computation reduces to polynomial root finding.

Explicit formulas for spectral measure derived.

Spectrum of compact perturbations is computable with error control.

Abstract

We address the computational spectral theory of Jacobi operators that are compact perturbations of the free Jacobi operator via the asymptotic properties of a connection coefficient matrix. In particular, for finite-rank perturbation we show that the computation of the spectrum can be reduced to a polynomial root finding problem, from a polynomial that is derived explicitly from the entries of a connection coefficient matrix. A formula for the spectral measure of the operator is also derived explicitly from these entries. The analysis is extended to trace-class perturbations. We address issues of computability in the framework of the Solvability Complexity Index, proving that the spectrum of compact perturbations of the free Jacobi operator is computable in finite time with guaranteed error control in the Hausdorff metric on sets.

Equations

$$J=\left(\begin{array}[]{cccc}\alpha_{0}&\beta_{0}&&\\ \beta_{0}&\alpha_{1}&\beta_{1}&\\ &\beta_{1}&\alpha_{2}&\smash{\raisebox{-4.0pt}{\footnotesize{$\ddots$}}}\\ &&\smash{\raisebox{-4.0pt}{\footnotesize{$\ddots$}}}&\smash{\raisebox{-4.0pt}{\footnotesize{$\ddots$}}}\end{array}\right),$$

$$UJ{U}^{\ast }[f](s)=sf(s),$$

$$P_k(s)=c_{0k}{Q}_0(s)+c_{1k}{Q}_1(s)+\cdots +c_{kk}{Q}_k(s).$$

$$\Delta=\left(\begin{array}[]{cccc}0&\frac{1}{2}&&\\ \frac{1}{2}&0&\frac{1}{2}&\\ &\frac{1}{2}&0&\smash{\raisebox{-4.0pt}{\footnotesize{$\ddots$}}}\\ &&\smash{\raisebox{-4.0pt}{\footnotesize{$\ddots$}}}&\smash{\raisebox{-4.0pt}{\footnotesize{$\ddots$}}}\end{array}\right),$$

$${\alpha }_k=0,{\beta }_{k-1}=\frac{1}{2}\text{ for all }k\ge n,$$

$$\sigma (J)=[-1,1]\cup \left\{\lambda (z_k):c(z_k)=0,z_k\in \mathbb{D}\right\},$$

$$\mathrm{d}\mu (s)=\frac{\mathrm{d}{\mu }_{\Delta }(s)}{|c(e^{i\theta }){|}^2}+\sum _{k=1}^r\frac{(z_k-z_k^{-1}{)}^2}{z_k{c}^{\prime }(z_k)c(z_k^{-1})}\mathrm{d}{\delta }_{\lambda (z_k)}(s),$$

$$\sum _{k=0}^{\infty }\left|{\alpha }_k\right|+\left|{\beta }_k-\frac{1}{2}\right|<\infty ,$$

$$CJ=DC,$$

$$\frac{\mathrm{d}\nu }{\mathrm{d}\mu }=\sum _{k=0}^{\infty }{c}_{0,k}{P}_k.$$

$$\Vert C{\Vert }_2^2=\underset{s\in \sigma (J)}{\mu -ess\sup }\left|\frac{\mathrm{d}\nu }{\mathrm{d}\mu }(s)\right|.$$

$$u_0(z)=\underset{k\to \infty }{\lim }(1-z^2)z^k{P}_k\left(\frac{1}{2}\left(z+z^{-1}\right)\right),z\in \mathbb{D}.$$

$$G(\lambda )=⟨e_0,(J-\lambda {)}^{-1}{e}_0⟩$$

$$⟨x,y⟩:=\sum _{k=0}^{\infty }{\bar{x}}_k{y}_k.$$

$$G(\lambda )=\int (s-\lambda {)}^{-1}\mathrm{d}\mu (s).$$

$$\frac{1}{2}\mu (\{s_1\})+\mu ((s_1,s_2))+\frac{1}{2}\mu (\{s_2\})=\underset{ϵ\searrow 0}{\lim }\frac{1}{\pi }{\int }_{s_1}^{s_2}\mathrm{Im}G(s+iϵ)\mathrm{d}s.$$

$$\sigma (J)=\mathrm{supp}(\mu )=\overline{\{s\in \mathbb{R}:\underset{ϵ\searrow 0}{\mathrm{liminf}}\mathrm{Im}G(s+iϵ)>0\}}.$$

$$\mu (\{s\})=\underset{ϵ\searrow 0}{\lim }\frac{ϵ}{i}G(s+iϵ)$$

$$\underset{ϵ\searrow 0}{\mathrm{liminf}}\mathrm{Im}G(s+iϵ)>0.$$

$$s{P}_k(s)$$

$$P_{-1}(s)$$

$$\int P_j(s)P_k(s)\mathrm{d}\mu (s)={\delta }_{jk}.$$

$$UJ{U}^{\ast }[f](s)=sf(s).$$

$$⟨e_i,f(J)e_j⟩=\int f(s)P_i(s)P_j(s)\mathrm{d}\mu (s).$$

$$\lambda P_k^{\mu }(\lambda )$$

$$P_0^{\mu }(\lambda )$$

$$P_k^{\mu }(\lambda )=\int \frac{P_k(s)-P_k(\lambda )}{s-\lambda }\mathrm{d}\mu (s),\lambda \in \mathbb{C}\setminus \sigma (J).$$

$$f^{\mu }(\lambda )=\int \frac{f(s)-f(\lambda )}{s-\lambda }\mathrm{d}\mu (s).$$

$$D=\left(\begin{array}[]{cccc}\gamma_{0}&\delta_{0}&&\\ \delta_{0}&\gamma_{1}&\delta_{1}&\\ &\delta_{1}&\gamma_{2}&\smash{\raisebox{-4.0pt}{\footnotesize{$\ddots$}}}\\ &&\smash{\raisebox{-4.0pt}{\footnotesize{$\ddots$}}}&\smash{\raisebox{-4.0pt}{\footnotesize{$\ddots$}}}\end{array}\right),$$

$$P_k(s)=c_{0k}{Q}_0(s)+c_{1k}{Q}_1(s)+\cdots +c_{kk}{Q}_k(s).$$

$$\left(\begin{array}[]{c}P_{0}(s)\\ P_{1}(s)\\ P_{2}(s)\\ \vdots\end{array}\right)=C^{T}\left(\begin{array}[]{c}Q_{0}(s)\\ Q_{1}(s)\\ Q_{2}(s)\\ \vdots\end{array}\right)\text{ for all }s\in\mathbb{C}.$$

Full text

Keywords   Jacobi operator, spectral measure, Solvability Complexity Index, orthogonal polynomials

AMS classification numbers   Primary 47B36, Secondary 47A10, 47A75, 65J10

1 Introduction

A Jacobi operator is a selfadjoint operator on $\ell^{2}=\ell^{2}(\{0,1,2,\ldots\})$, which with respect to the standard orthonormal basis $\{e_{0},e_{1},e_{2},\ldots\}$ has a tridiagonal matrix representation,

[TABLE]

where $\alpha_{k}$ and $\beta_{k}$ are real numbers with $\beta_{k}>0$. In the special case where $\{\beta_{k}\}_{k=0}^{\infty}$ is a constant sequence, this operator can be interpreted as a discrete Schrödinger operator on the half line, and its spectral theory is arguably its most important aspect.

The spectral theorem for Jacobi operators guarantees the existence of a probability measure $\mu$ supported on the spectrum $\sigma(J)\subset\mathbb{R}$, called the spectral measure, and a unitary operator $U:\ell^{2}\to L^{2}_{\mu}(\mathbb{R})$ such that

[TABLE]

for all $f\in L^{2}_{\mu}(\mathbb{R})$ [16]. The coefficients $\{\alpha_{k}\}_{k\geq 0}$ and $\{\beta_{k}\}_{k\geq 0}$ are the three–term recurrence coefficients of the orthonormal polynomials $\{P_{k}\}_{k\geq 0}$ with respect to $\mu$, which are given by $P_{k}=Ue_{k}$.

Suppose we have a second Jacobi operator, $D$, for which the spectral theory is known analytically. The point of this paper is to show that for certain classes of Jacobi operators $J$ and $D$, the computation and theoretical study of the spectrum and spectral measure of $J$ can be conducted effectively using the connection coefficient matrix between $J$ and $D$, combined with known properties of $D$.

Definition 1.1**.**

The connection coefficient matrix $C=C_{J\to D}=(c_{ij})_{i,j=0}^{\infty}$ is the upper triangular matrix representing the change of basis between the orthonormal polynomials $(P_{k})_{k=0}^{\infty}$ of $J$, and the orthonormal polynomials $(Q_{k})_{k=0}^{\infty}$ of $D$, whose entries satisfy,

[TABLE]

We pay particular attention to the case where $D$ is the so-called free Jacobi operator,

[TABLE]

and $J$ is a Jacobi operator of the form $J=\Delta+K$, where $K$ is compact. In the discrete Schrödinger operator setting, $K$ is a diagonal potential function which decays to zero at infinity [41]. Another reason this class of operators is well studied is because the Jacobi operators for the classical Jacobi polynomials are of this form [36]. Since scaling and shifting by the identity operator affects the spectrum in a trivial way, results about these Jacobi operators $J$ apply to all Jacobi operators which are compact perturbations of a Toeplitz operator (Toeplitz-plus-compact).

The spectral theory of Toeplitz operators such as $\Delta$ is known explicitly [7]. The spectral measure of $\Delta$ is the semi-circle $\mathrm{d}\mu_{\Delta}(s)=\frac{2}{\pi}(1-s^{2})^{\frac{1}{2}}$ (restricted to $[-1,1]$), and the orthonormal polynomials are the Chebyshev polynomials of the second kind, $U_{k}(s)$. We prove the following new results.

If $J$ is a Jacobi operator which is a finite rank perturbation of $\Delta$ (Toeplitz-plus-finite-rank), i.e. there exists $n$ such that

[TABLE]

• •

Theorem 4.8: The connection coefficient matrix $C_{J\to\Delta}$ can be decomposed into $C_{\rm Toe}+C_{\rm fin}$ where $C_{\rm Toe}$ is Toeplitz, upper triangular and has bandwidth $2n-1$, and the entries of $C_{\rm fin}$ are zero outside the $n-1\times 2n-1$ principal submatrix.

• •

Theorem 4.22: Let $c$ be the Toeplitz symbol of $C_{\rm Toe}$. It is a degree $2n-1$ polynomial with $r\leq n$ roots inside the complex open unit disc $\mathbb{D}$, all of which are real and simple. The spectrum of $J$ is

[TABLE]

where $\lambda(z):=\frac{1}{2}(z+z^{-1}):\mathbb{D}\to\mathbb{C}\setminus[-1,1]$ is the Joukowski transformation. The spectral measure of $J$ is given by the formula

[TABLE]

where $\cos(\theta)=s$. The denominator in the first term can be expressed using the polynomial, $|c(e^{i\theta})|^{2}=\sum_{k=0}^{2n-1}\langle C^{T}e_{k},C^{T}e_{0}\rangle U_{k}(s)$.

For $R>0$, define the Banach space $\ell^{1}_{R}$ to be scalar sequences such that $\sum_{k=0}^{\infty}|v_{k}|R^{k}<\infty$. If $J$ is a trace class perturbation of $\Delta$ (Toeplitz-plus-trace-class), i.e.,

[TABLE]

• •

Theorem 5.11: $C=C_{J\to\Delta}$ is bounded as an operator from $\ell^{1}_{R}$ into itself, for all $R>1$. Further, we have the decomposition $C=C_{\rm Toe}+C_{K}$ where $C_{\rm Toe}$ is upper triangular Toeplitz and $C_{K}$ is compact as an operator from $\ell^{1}_{R}$ into itself for all $R>1$.

• •

Theorem 5.13 and Theorem 5.15 : The Toeplitz symbol of $C_{\rm Toe}$, $c$, is analytic in the complex unit disc. The discrete eigenvalues, as in the Toeplitz-plus-finite-rank case are of the form $\frac{1}{2}(z_{k}+z_{k}^{-1})$ where $z_{k}$ are the roots of $c$ in the open unit disc.

The relevance of the space $\ell^{1}_{R}$ here is that for an upper triangular Toeplitz matrix which is bounded as an operator from $\ell^{1}_{R}$ to itself for all $R>1$, the symbol of that operator is analytic in the open unit disc.

Following the pioneering work of Ben-Artzi–Colbrook–Hansen–Nevanlinna–Seidel on the Solvability Complexity Index [4, 5, 28, 9, 8], we prove two theorems about computability. We assume real number arithmetic, and the results do not necessarily apply to algorithms using floating point arithmetic.

• •

Theorem 6.8: If $J$ is a Toeplitz-plus-finite-rank Jacobi operator, then in a finite number of operations, the absolutely continuous part of the spectral measure is computable exactly, and the locations and weights of the discrete part of the spectral measure are computable to any desired accuracy. If the rank of the perturbation is known a priori then the algorithm can be designed to terminate with guaranteed error control.

• •

Theorem 6.10: If $J=\Delta+K$ is a Toeplitz-plus-compact Jacobi operator, then in a finite number of operations, the spectrum of $J$ is computable to any desired accuracy in the Hausdorff metric on subsets of $\mathbb{R}$. If the quantity $\sup_{k\geq m}|\alpha_{k}|+\sup_{k\geq m}|\beta_{k}-\frac{1}{2}|$ can be estimated for all $m$, then the algorithm can be designed to terminate with guaranteed error control.

The present authors consider these results to be the beginning of a long term project on the computation of spectra of structured operators. Directions for future research are outlined in Section 7.

1.1 Relation to existing work on connection coefficients

Recall that $J$ and $D$ are Jacobi operators with orthonormal polynomials $\{P_{k}\}_{k=0}^{\infty}$ and $\{Q_{k}\}_{k=0}^{\infty}$ respectively, and spectral measures $\mu$ and $\nu$ respectively. Uvarov gave expressions for the relation between $\{P_{k}\}_{k=0}^{\infty}$ and $\{Q_{k}\}_{k=0}^{\infty}$ in the case where the Radon–Nikodym derivative $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$ is a rational function, utilising connection coefficients [47, 48]. Decades later, Kautsky and Golub related properties of $J$, $D$ and the connection coefficients matrix $C=C_{J\to D}$ with the Radon-Nikodym derivative $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$, with applications to practical computation of Gauss quadrature rules in mind [32]. In the setting of Gauss quadrature rules, the connection coefficients are usually known as modified or mixed moments, (see [23, 40, 55]). The following results, which we prove in Section 3 for completeness, are straightforward generalisations of what can be found in the papers cited in this paragraph from the 1960’s, 1970’s and 1980’s.

The Jacobi operators and the connection coefficients satisfy

[TABLE]

which makes sense as operators acting on finitely supported sequences (this is made clear and proved in Theorem 3.4). A finite-dimensional version of this result with a rank-one remainder term first appears in [32, Lem. 1]. The connection coefficients matrix also determines the existence and certain properties of the Radon–Nikodym derivative $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$:

• •

Proposition 3.6: $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}\in L^{2}_{\mu}(\mathbb{R})$ if and only if the first row of $C$ is an $\ell^{2}$ sequence, in which case

[TABLE]

• •

Proposition 3.9: If $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}\in L^{\infty}_{\mu}(\mathbb{R})$ then $C$ is a bounded operator on $\ell^{2}$ and

[TABLE]

• •

Corollary 3.10: If both $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}\in L^{\infty}_{\mu}(\mathbb{R})$ and $\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\in L^{\infty}_{\nu}(\mathbb{R})$ then $C$ is bounded and invertible on $\ell^{2}$.

1.2 Relation to existing work on spectra of Jacobi operators

There has been extensive work in the last 20 years or so on the spectral theory of Jacobi operators which are compact perturbations of the free Jacobi operator, particularly with applications to quantum theory and random matrix theory in mind [33, 14, 15, 42, 20, 22]. The literature focuses on so-called sum rules, which are certain remarkable relationships between the spectral measures and the entries of Jacobi operators, and builds upon some 20th century developments in the Szegő theory of orthogonal polynomials [44, 30, 25, 18, 52, 49, 50, 35, 51].

For a Jacobi operator $J$, which is a compact perturbation of the free Jacobi operator $\Delta$, with orthogonal polynomials $\{P_{k}\}_{k=0}^{\infty}$ and spectral measure $\mathrm{d}\mu$, the central analytical object of interest is the function

[TABLE]

Conditions on the perturbation $J-\Delta$ yield a nicer $u_{0}$. For instance, when $J-\Delta$ is trace class, $u_{0}$ is analytic and its roots in the unit disc are in one-to-one correspondence with the eigenvalues of $J$, and when $J-\Delta$ is of finite rank, $u_{0}$ is a polynomial [33].

The function $u_{0}$ defined as in (1.12) can be arrived at through several different definitions: the Jost functions, the perturbation determinant, the Szegő function and the Geronimo–Case functions (see the introduction of [33]). One contribution of the present paper is that it provides a new interpretation of $u_{0}$: the Toeplitz symbol $c(z)$. Furthermore, this new interpretation also has a matrix associated to it, the connection coefficients matrix $C$, which is defined for any two Jacobi operators, not only perturbations of the free Jacobi operator. This opens the door to generalisation, something the present authors intend to pursue in the future.

Very recently, Colbrook, Bogdan, and Hansen have introduced techniques for computing spectra with error control [10] which work on quite broad classes of operators. Colbrook [11] extended these techniques to computing spectral measures of operators, including Jacobi operators as a special case. While this work applies to broader classes of operators than ours, it gives less accurate results for the class of Jacobi operators we consider, and in particular does not produce explicit formulae such as Theorem 4.22 or as precise results as Theorem 6.10. In particular, our assumptions on the structure of the operator are sufficient to produce better results in terms the SCI hierarchy: we can compute the discrete spectra and spectral measures of trace-class perturbations or compact perturbations with known decay with error control ($\Delta_{1}$ in the notation of [11]), the spectral measure of Jacobi operators that are finite rank perturbations of $\Delta$ in finite operations ($\Delta_{0}$), and the absolutely continuous spectrum is always $[-1,1]$. This compares favourably to [11] which proves $\Delta_{2}$ classification results (one limit with no error control) for the spectral measures, projections, functional calculus and Radon-Nikodym derivatives of a larger class of operators. Furthermore, the generality of [11] means that the location of the absolutely continuous and pure point spectra are no longer known, as is the case for our class of operators. This causes the computation of spectral decompositions to become very difficult and higher up in the SCI hierarchy.

1.3 Outline of the paper

• •

In Section 2 we outline basic, established results about spectral theory of Jacobi operators.

• •

In Section 3 we discuss the basic properties of the connection coefficients matrix $C_{J\to D}$ for general Jacobi operators $J$ and $D$, and how they relate to the spectra of $J$ and $D$.

• •

In Section 4 we show how connection coefficient matrices apply to Toeplitz-plus-finite-rank Jacobi operators, and in Section 5 we extend these results to the Toeplitz-plus-trace-class case.

• •

Section 6 is devoted to issues of computability.

• •

Appendix A gives an array of numerical examples produced using an open source Julia package SpectralMeasures.jl [53] that implements the ideas of this paper. It makes extensive use of the open source Julia package ApproxFun.jl [37, 38], in particular the features for defining and manipulating functions and infinite-dimensional operators.

Acknowledgements

The authors thank the following people for helpful discussions on this work as it developed: Anders Hansen, Arieh Iserles, Andrew Swan, Peter Clarkson, Alfredo Deaño, Alex Townsend, Walter Van Assche, Bernhard Beckermann, Matthew Colbrook, Giles Shaw and David Sanders, amongst others. We are also grateful to the anonymous referees for their comments which greatly improved the quality of the paper. The first author was supported by the London Mathematical Society Cecil King Travel Scholarship to visit the University of Sydney in 2016, the UK Engineering and Physical Sciences Research Council (EPSRC) grant EP/H023348/1 for the University of Cambridge Centre for Doctoral Training, the Cambridge Centre for Analysis, by the FWO research project G.A004.14 at KU Leuven, Belgium, and by an FWO Postdoctoral Fellowship. The second author was supported by a Leverhulme Trust Research Project Grant. Portions of this manuscript were submitted in an earlier form as part of the PhD thesis of first author at the University of Cambridge in 2017.

2 Spectral theory of Jacobi operators

In this section we present well known results about the spectra of Jacobi operators. This gives a self-contained account of what is required to prove the results later in the paper, and sets the notation.

Definition 2.1**.**

Define the principal resolvent function for $\lambda\in\mathbb{C}\setminus\sigma(J)$,

[TABLE]

where

[TABLE]

Theorem 2.2** ([16, 45, 41]).**

Let $J$ be a bounded Jacobi operator.

1. (i)

There exists a unique compactly supported probability measure $\mu$ on $\mathbb{R}$, called the spectral measure of $J$, such that

[TABLE] 2. (ii)

For any $s_{1}<s_{2}$ in $\mathbb{R}$,

[TABLE] 3. (iii)

The spectrum of $J$ is

[TABLE]

The point spectrum $\sigma_{p}(J)$ of $J$ is the set of points $s\in\mathbb{R}$ such that the limit

[TABLE]

exists and is positive.

The continuous spectrum of $J$ is the set of points $s\in\mathbb{R}$ such that $\mu(\{s\})=0$ but

[TABLE]

The measure $\mu$ is the spectral measure that appears in the spectral theorem for self-adjoint operators on Hilbert space [21], as demonstrated by the following theorem [16, 45, 41].

Definition 2.3**.**

The orthonormal polynomials for $J$ are $P_{0},P_{1},P_{2},\ldots$ defined by the three term recurrence

[TABLE]

Theorem 2.4** ([16]).**

Let $J$ be a bounded Jacobi operator and let $P_{0},P_{1},P_{2},\ldots$ be as defined in Definition 2.3. Then we have the following.

1. (i)

The polynomials are such that $P_{k}(J)e_{0}=e_{k}$. 2. (ii)

The polynomials are orthonormal with respect to the spectral measure of $J$,

[TABLE] 3. (iii)

Define the unitary operator $U:\ell^{2}\to L^{2}_{\mu}(\mathbb{R})$ such that $Ue_{k}=P_{k}$. Then for all $f\in L^{2}_{\mu}(\mathbb{R})$,

[TABLE] 4. (iv)

For all polynomials $f$, the entries of $f(J)$ are equal to,

[TABLE]

For $f\in L^{1}_{\mu}(\mathbb{R})$, this formula defines the matrix $f(J)$.

The following definition is standard in orthogonal polynomial theory.

Definition 2.5** ([24, 50]).**

The first associated polynomials for $J$ are $P^{\mu}_{0},P^{\mu}_{1},P^{\mu}_{2},\ldots$ defined by the three term recurrence

[TABLE]

The relevance of the first associated polynomials for this work is the following integral formula.

Lemma 2.6**.**

([24, pp. 17,18]) The first associated polynomials are given by the integral formula

[TABLE]

For notational convenience we also define the $\mu$-derivative of a general polynomial.

Definition 2.7**.**

Let $\mu$ be a probability measure compactly supported on the real line and let $f$ be a polynomial. The $\mu$-derivative of $f$ is the polynomial defined by

[TABLE]

3 Connection coefficient matrices

In this section we give preliminary results to indicate the relevance of connection coefficient matrices to spectral theory of Jacobi operators.

3.1 Basic properties

As in the introduction, consider a second bounded Jacobi operator,

[TABLE]

with principal resolvent function $H(z)$, spectral measure $\nu$ and orthogonal polynomials denoted $Q_{0},Q_{1},Q_{2},\ldots$. In the introduction (Definition 1.1) we defined the connection coefficient matrix between $J$ and $D$, $C=C_{J\to D}$ to have entries satisfying

[TABLE]

Definition 3.1**.**

Denote the space of complex-valued sequences with finitely many nonzero elements by $\ell_{\mathcal{F}}$, and its algebraic dual, the space of all complex-valued sequences, by $\ell_{\mathcal{F}}^{\star}$.

Note that $C:\ell_{\mathcal{F}}\rightarrow\ell_{\mathcal{F}}$, because it is upper triangular, and $C^{T}:\ell_{\mathcal{F}}^{\star}\rightarrow\ell_{\mathcal{F}}^{\star}$, because it is lower triangular, and thus we may write

[TABLE]

By orthonormality of the polynomial sequences the entries can also be interpreted as

[TABLE]

where $\langle\cdot,\cdot\rangle_{\nu}$ is the standard inner product on $L^{2}_{\nu}(\mathbb{R})$.

A recurrence relationship for the connection coefficients matrix was discovered by Sack and Donovan [40] and independently by Wheeler [55], in the context of Gauss quadrature formulae.

Lemma 3.2**.**

[40, 55]** The entries of the connection coefficients matrix $C_{J\to D}$ satisfy the following 5-point discrete system:

[TABLE]

*with boundary conditions *

[TABLE]

Proof.

Assume by convention that $c_{ij}=0$ if $i=-1$ or $j=-1$. Now using this boundary condition and the three term recurrences for the polynomial sequences, we see that

[TABLE]

and

[TABLE]

Since $\langle sQ_{i}(s),P_{j}(s)\rangle_{\nu}=\langle Q_{i}(s),sP_{j}(s)\rangle_{\nu}$, we have the result for the interior points $0\leq i<j$.

The remaining boundary conditions come from $c_{i0}=\langle Q_{i},P_{0}\rangle_{\nu}$ which equals 1 if $i=0$ and 0 otherwise. ∎

The 5-point recurrence formula can be restated as infinite-vector-valued three-term recurrence relations for rows and columns of $C$.

Corollary 3.3**.**

The columns of $C$ satisfy

[TABLE]

Consequently the $j$th column can be written $c_{*,j}=P_{j}(D)e_{0}$.

The rows of $C$ satisfy

[TABLE]

Consequently, the $i$th row can be written $c_{i,*}=c_{0,*}Q_{i}(J)$.

Proof.

The 5-point discrete system described in Lemma 3.2 can be used to find an explicit linear recurrence to compute the entries of $C$,

[TABLE]

The recurrences of the rows and columns of $C$ are these written in vectorial form.

The consequences follow from the uniqueness of solution to second order difference equations with two initial data (adding $c_{-1,*}=0$ and $c_{*,-1}=0$). ∎

3.2 Connection coefficients and spectral theory

The following theorems give precise results about how the connection coefficients matrix $C$ can be useful for studying and computing the spectra of Jacobi operators.

Theorem 3.4** ([32]).**

Let $J$ and $D$ be bounded Jacobi operators and $C=C_{J\to D}$ the connection coefficients matrix. For all polynomials $p$, we have the following as operators from $\ell_{\mathcal{F}}$ to $\ell_{\mathcal{F}}$,

[TABLE]

Remark 3.5*.*

This is a generalisation of the result of Kausky and Golub [32, Lem 1], that

[TABLE]

where $C_{N\times N},J_{N\times N},D_{N\times N}$ are the principal $N\times N$ submatrices of $C,J,D$, and $c_{N}$ is a certain vector in $\mathbb{R}^{N}$.

Proof.

First we begin with the case $p(z)=z$. By definition,

[TABLE]

Then by Corollary 3.3, this is equal to $Dc_{*,0}$, which is equal to $DCe_{0}$. Now, for any $j>0$,

[TABLE]

Then by Corollary 3.3, this is equal to $Dc_{*,j}$, which is equal to $DCe_{j}$. Hence $CJ=DC$.

Now, when $f(z)=z^{k}$ for any $k>0$, $D^{k}C=D^{k-1}CJ=\cdots=CJ^{k}$. By linearity $Cf(J)=f(D)C$ for all polynomials $f$. ∎

We believe that the basic properties relating the Radon-Nikodym derivative $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$ and the connection coefficients matrix $C$ have not been given in the literature before, but follow naturally from discussions in, for example, [27, Ch. 5], on mixed moments and modifications of weight functions for orthogonal polynomials.

Proposition 3.6**.**

Let $J$ and $D$ be bounded Jacobi operators with spectral measures $\mu$ and $\nu$ respectively, and connection coefficient matrix $C=C_{J\to D}$. Then

[TABLE]

in which case

[TABLE]

Proof.

Suppose first that $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}\in L^{2}_{\mu}(\mathbb{R})$. Then $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}=\sum_{k=0}^{\infty}a_{k}P_{k}$, for some $a\in\ell^{2}$, because $P_{0},P_{1},P_{2},\ldots$ is an orthonormal basis of $L^{2}_{\mu}(\mathbb{R})$. The coefficients are given by,

[TABLE]

Hence $c_{0,*}\in\ell^{2}$ and gives the $P_{k}$ coefficients of $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$.

Conversely, suppose that $c_{0,*}\in\ell^{2}$. Then the function $\sum_{k=0}^{\infty}c_{0,k}P_{k}$ is in $L^{2}_{\mu}(\mathbb{R})$, and by the same manipulations as above its projections onto polynomial subspaces are equal to that of $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$. ∎

Remark 3.7*.*

If we have a situation in which $c_{0,*}\in\ell^{2}$, we can by Proposition 3.6 and the existence of the Radon–Nikodym derivative on $\mathrm{supp}(\mu)$ deduce that $\sigma(D)\subset\sigma(J)$ and the function defined by $\sum_{k=0}^{\infty}c_{0,k}P_{k}$ is zero on $\sigma(J)\setminus\sigma(D)$. This observation translates into a rootfinding problem in Section 4.

Lemma 3.8**.**

Let $J$ and $D$ be bounded Jacobi operators with spectral measures $\mu$ and $\nu$ respectively, and connection coefficient matrix $C=C_{J\to D}$. If $\nu$ is absolutely continuous with respect to $\mu$, then as operators mapping $\ell_{\mathcal{F}}\to\ell_{\mathcal{F}}^{\star}$,

[TABLE]

Here the matrix $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}(J)$ is interpreted as in Theorem 2.4 part (iii).

Proof.

Note first that since $C:\ell_{\mathcal{F}}\to\ell_{\mathcal{F}}$ and $C^{T}:\ell_{\mathcal{F}}^{\star}\to\ell_{\mathcal{F}}^{\star}$, $C^{T}C$ is well-defined $\ell_{\mathcal{F}}\to\ell_{\mathcal{F}}^{\star}$. Then we have,

[TABLE]

This completes the proof. ∎

Proposition 3.9**.**

Let $J$ and $D$ be bounded Jacobi operators with spectral measures $\mu$ and $\nu$ respectively, and connection coefficient matrix $C=C_{J\to D}$. If $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}\in L^{\infty}_{\mu}(\mathbb{R})$ then $C$ is a bounded operator on $\ell^{2}$ and

[TABLE]

Here $\|\cdot\|_{2}$ is the operator norm from $\ell^{2}\to\ell^{2}$ and $\operatorname*{\mu-ess\,sup}$ is the supremum up to $\mu$-almost everywhere equivalence of functions.

Proof.

Since $\mu$ is a probability measure (and hence $\sigma$-finite), we have the standard characterisation,

[TABLE]

which can be modified to

[TABLE]

by associating positive functions $g\in L^{1}_{\mu}(\mathbb{R})$ with their square-roots $f\in L_{\mu}^{2}(\mathbb{R})$. Now, since $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}\in L_{\mu}^{\infty}(\mathbb{R})$, this supremum can actually be taken over all polynomials by the following argument. If $f\in L_{\mu}^{2}(\mathbb{R})$ with $\|f\|_{2,\mu}\leq 1$ then for any $\epsilon>0$ there exists polynomial $p$ such that $\|p\|_{2,\mu}\leq 1$ and $\|f-p\|_{2}\leq\epsilon$, since polynomials are dense in $L_{\mu}^{2}(\mathbb{R})$ (this follows, for example, from the compact support of $\mu$ [1]). It is readily shown using Hölder and triangle inequalities that

[TABLE]

so that any supremum can be arbitrarily approximated using a polynomial.

Since $\{P_{k}\}_{k=0}^{\infty}$ is a complete orthonormal basis of $L_{\mu}^{2}(\mathbb{R})$ (completeness holds as a result of $\mu$ having compact support), then $f\in L_{\mu}^{2}(\mathbb{R})$ if and only if there is a unique sequence $v\in\ell^{2}$ such that $f=\sum_{k=0}^{\infty}v_{k}P_{k}$ with the series converging in the $L_{\mu}^{2}(\mathbb{R})$ norm. Furthermore, $\|f\|_{2}=\|v\|_{2}$. Hence,

[TABLE]

Since $\ell_{\mathcal{F}}$ is a dense subspace of $\ell^{2}$, we have $\|C\|_{2}=\sup_{\begin{subarray}{c}v\in\ell_{\mathcal{F}}\end{subarray},\\ \|v\|_{2}=1}\|Cv\|_{2}$. Now, $\|Cv\|_{2}^{2}=\langle v,C^{T}Cv\rangle$, and by Lemma 3.8, $C^{T}C=\frac{\mathrm{d}\nu}{\mathrm{d}\mu}(J)$. Therefore,

[TABLE]

By Theorem 2.4 part (iv), we conclude that

[TABLE]

Therefore, $\|C\|_{2}^{2}=\operatorname*{\mu-ess\,sup}_{s\in\sigma(J)}\left|\frac{\mathrm{d}\nu}{\mathrm{d}\mu}(s)\right|$ and $C$ is bounded. ∎

Corollary 3.10**.**

Let $J$ and $D$ be bounded Jacobi operators with spectral measures $\mu$ and $\nu$ respectively, and connection coefficient matrix $C=C_{J\to D}$. If $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}\in L^{\infty}_{\mu}(\mathbb{R})$ and $\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\in L^{\infty}_{\nu}(\mathbb{R})$ then $C$ is bounded and invertible on $\ell^{2}$.

Proof.

By Proposition 3.9, $C_{J\to D}$ is bounded if $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}\in L^{\infty}_{\mu}(\mathbb{R})$, and $C_{D\to J}$ is bounded if $\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\in L^{\infty}_{\nu}(\mathbb{R})$. Combining this the fact that $C_{J\to D}^{-1}=C_{D\to J}$, as operators from $\ell_{\mathcal{F}}$ to itself, we complete the proof. ∎

The following definition and lemma are useful later.

Definition 3.11**.**

Given polynomial sequences $P_{0},P_{1},P_{2},\ldots$ and $Q_{0},Q_{1},Q_{2},\ldots$ for Jacobi operators $J$ and $D$ respectively, we define the matrix $C^{\mu}$ to be the connection coefficients matrix between $P_{0}^{\mu},P_{1}^{\mu},P_{2}^{\mu},\ldots$ and $Q_{0},Q_{1},Q_{2},\ldots$ as in Definition 1.1, where $P_{0}^{\mu},P_{1}^{\mu},P_{2}^{\mu},\ldots$ are the first associated polynomials for $J$ as in Definition 2.5. Noting that the lower triangular matrix $(C^{\mu})^{T}$ is a well defined operator from $\ell_{\mathcal{F}}^{\star}$ into itself, we have

[TABLE]

Remark 3.12*.*

Note that $C^{\mu}$ is strictly upper triangular, because the first associated polynomials have their degrees one less than their indices.

Lemma 3.13**.**

The operator $C^{\mu}$ as defined above for $C_{J\to D}$ is in fact $\beta_{0}^{-1}\left(0,C_{J^{\mu}\to D}\right)$, where

[TABLE]

Proof.

The (unique) orthonormal polynomials for $J^{\mu}$ are $\beta_{0}P_{1}^{\mu},\beta_{0}P_{2}^{\mu},\beta_{0}P_{3}^{\mu},\ldots$, and $P^{\mu}_{0}=0$. ∎

4 Toeplitz-plus-finite-rank Jacobi operators

In this section we present several novel results which show how the connection coefficient matrices can be used for computing the spectral measure of a Toeplitz-plus-finite-rank Jacobi operator.

4.1 Jacobi operators for Chebyshev polynomials

There are two particular Jacobi operators with Toeplitz-plus-finite-rank structure that are of great interest,

[TABLE]

The spectral measures of $\Delta$ and $\Gamma$ are

[TABLE]

supported on $[-1,1]$.

Using results of Stieltjes in his seminal paper [43], [1, App.], the principal resolvent can be written elegantly as a continued fraction,

[TABLE]

Using this gives explicit expressions for the principal resolvents,

[TABLE]

Remark 4.1*.*

We must be careful about which branch we refer to when we write the resolvents in this explicit form. Wherever $\sqrt{}$ is written above we mean the standard branch that is positive on $(0,\infty)$ with branch cut $(-\infty,0]$. This gives a branch cut along $[-1,1]$ in both cases, the discontinuity of $G$ across which makes the Perron–Stieltjes inversion formula in Theorem 2.2 work. It also ensures the $\mathcal{O}(\lambda^{-1})$ decay resolvents enjoy as $\lambda\to\infty$.

The orthonormal polynomials for $\Delta$ are the Chebyshev polynomials of the second kind, which we denote $U_{k}(s)$,

[TABLE]

The orthonormal polynomials for $\Gamma$ are the normalised Chebyshev polynomials of the first kind, which we denote $\tilde{T}_{k}(s)$. Note that these are not the usual Chebyshev polynomials of the first kind (denoted $T_{k}(s)$) [24, 16]. We in fact have,

[TABLE]

The first associated polynomials have simple relationships with the orthonormal polynomials,

[TABLE]

4.2 Basic perturbations

In this section we demonstrate for two simple, rank-one perturbations of $\Delta$ how the connection coefficient matrix relates properties of the spectrum of the operators. This will give some intuition as to what to expect in more general cases.

Example 4.2** (Basic perturbation 1).**

Let $\alpha\in\mathbb{R}$, and define

[TABLE]

Then the connection coefficient matrix $C_{J_{\alpha}\to\Delta}$ is the bidiagonal Toeplitz matrix

[TABLE]

This can be computed using the explicit recurrences (3.3)–(3.7). The connection coefficient matrix $C_{\Delta\to J_{\alpha}}$ (which is the inverse of $C_{J_{\alpha}\to\Delta}$ on $\ell_{\mathcal{F}}$) is the full Toeplitz matrix

[TABLE]

From this we see that $C=C_{J_{\alpha}\to\Delta}$ has a bounded inverse in $\ell^{2}$ if and only if $|\alpha|<1$. Hence by Theorem 3.4, if $|\alpha|<1$ then $CJ_{\alpha}C^{-1}=\Delta$ with each operator bounded on $\ell^{2}$, so that $\sigma(J_{\alpha})=\sigma(\Delta)=[-1,1]$. We will discuss what happens when $|\alpha|\geq 1$ later in the section.

Example 4.3** (Basic perturbation 2).**

Let $\beta>0$, and define

[TABLE]

Then the connection coefficient matrix $C_{J_{\beta}\to\Delta}$ is the banded Toeplitz-plus-rank-1 matrix

[TABLE]

Just as in Example 4.2, this can be computed using the explicit recurrences (3.3)–(3.7). The connection coefficient matrix $C_{\Delta\to J_{\beta}}$ (which is the inverse of $C_{J_{\beta}\to\Delta}$ on $\ell_{\mathcal{F}}$) is the Toeplitz-plus-rank-1 matrix

[TABLE]

From this we see that $C=C_{J_{\beta}\to\Delta}$ has a bounded inverse on $\ell^{2}$ if and only if $\beta<\sqrt{2}$. Hence by Theorem 3.4, if $\beta<\sqrt{2}$ then $CJ_{\beta}C^{-1}=\Delta$ with each operator bounded on $\ell^{2}$, so that $\sigma(J_{\beta})=\sigma(\Delta)=[-1,1]$. We will discuss what happens when $\beta\geq\sqrt{2}$ later in the section. Note that the case $\beta=\sqrt{2}$ gives the Jacobi operator $\Gamma$ in equation (4.1).

4.3 Fine properties of the connection coefficients

The two basic perturbations of $\Delta$ discussed above give connection coefficient matrices that are highly structured. The following lemmata and theorems prove that this is no coincidence; in fact, if Jacobi operator $J$ is a finite-rank perturbation of $\Delta$ then $C_{J\to\Delta}$ is also a finite-rank perturbation of Toeplitz.

Remark 4.4*.*

Note for the following results that all vectors and matrices are indexed starting from 0.

Lemma 4.5**.**

If $\delta_{j}=\beta_{j}$ for $j\geq n$ then $c_{jj}=c_{nn}$ for all $j\geq n$.

Proof.

By the recurrence in Lemma 3.2, $c_{jj}=(\delta_{j-1}/\beta_{j-1})c_{j-1,j-1}$. The result follows by induction. ∎

Lemma 4.6**.**

Let $J$ and $D$ be Jacobi operators with coefficients $\{\alpha_{k},\beta_{k}\}$ and $\{\gamma_{k},\delta_{k}\}$ respectively, such that there exists an $n$ such that111More intuitively, the entries of $J$ and $D$ are both equal and Toeplitz, except in the principal $n\times n$ submatrix, where neither statement necessarily holds.

[TABLE]

Then the entries of the connection coefficient matrix $C=C_{J\to D}$ satisfy

[TABLE]

Remark 4.7*.*

This means that $C$ is of the form $C=C_{\rm Toe}+C_{\rm fin}$ where $C_{\rm Toe}$ is Toeplitz and $C_{\rm fin}$ is zero except in the first $n-1$ rows. For example, when $n=4$, we have the following structure

[TABLE]

Proof.

We prove by induction on $k=0,1,2,\ldots$ that

[TABLE]

We use the recurrences in Lemma 3.2 and equations (3.3)–(3.7). The base case $k=0$ is proved in Lemma 4.5. Now we deal with the second base case, $k=1$. For any $i\geq n$, we have $\beta_{i}=\delta_{i}=\beta_{i-1}=\delta_{i-1}$, and $\alpha_{i}=\gamma_{i}$, so

[TABLE]

Now we deal with the case $k>1$. For any $i\geq n$, we have $\delta_{i}=\delta_{i-1}=\beta_{i+k-2}=\beta_{i+k-1}$, and $\alpha_{i+k-1}=\gamma_{i}$, so

[TABLE]

The last line follows from the induction hypothesis for the case $k-2$ (hence why we needed two base cases). ∎

The special case in which $D$ is Toeplitz gives even more structure to $C$, as demonstrated by the following theorem. We state the results for a finite-rank perturbation of the free Jacobi operator $\Delta$, but they apply to general Toeplitz-plus-finite rank Jacobi operators because the connection coefficients matrix $C$ is unaffected by a scaling and shift by the identity applied to both $J$ and $D$.

Theorem 4.8**.**

Let $J$ be a Jacobi operator such that there exists an $n$ such that

[TABLE]

i.e. it is equal to the free Jacobi operator $\Delta$ outside the $n\times n$ principal submatrix. Then the entries of the connection coefficient matrix $C=C_{J\to\Delta}$ satisfy

[TABLE]

Remark 4.9*.*

This means that $C$ is of the form $C=C_{\rm Toe}+C_{\rm fin}$ where $C_{\rm Toe}$ is Toeplitz with bandwidth $2n-1$ and $C_{\rm fin}$ zero except for entries in the $(n-1)\times(2n-2)$ principal submatrix. For example when $n=4$, we have the following structure,

[TABLE]

Proof.

First we prove (4.7). Fix $i,j$ such that $i+j\geq 2n$. Note that the case $i\geq n$ is proven in Lemma 4.6. Hence we assume $i<n$, and therefore $j>n$. Using Lemma 3.2 and equations (3.3)–(3.7) we find the following recurrence. Substituting $\delta_{i}=\frac{1}{2}$, $\gamma_{i}=0$ for all $i$, and $\alpha_{k}=0$, $\beta_{k-1}=\frac{1}{2}$ for $k\geq n$ into the recurrence, we have

[TABLE]

Repeating this process on $c_{i+1,j-1}$ in the above expression gives

[TABLE]

Repeating the process on $c_{i+2,j-2}$ and so on eventually gives

[TABLE]

By Lemma 4.6, $c_{n,i+j-n}=c_{n-1,i+j-n-1}$, so we are left with $c_{i,j}=c_{i-1,j-1}$. This completes the proof of (4.7).

Now we prove (4.8). Let $j\geq 2n$. Then

[TABLE]

This is equal to zero by (4.7), because $1+(j-1)\geq 2n$. ∎

Corollary 4.10**.**

Let $C^{\mu}$ be as defined in Definition 3.11 for $C$ as in Theorem 4.8. Then $C^{\mu}=C^{\mu}_{\rm Toe}+C^{\mu}_{\rm fin}$, where $C^{\mu}_{\rm Toe}$ is Toeplitz with bandwidth $2n-2$ and $C^{\mu}_{\rm fin}$ is zero outside the $(n-2)\times(2n-1)$ principal submatrix.

Proof.

This follows from Theorem 4.8 applied to $J^{\mu}$ as defined in Lemma 3.13. ∎

Remark 4.11*.*

A technical point worth noting for use in proofs later is that for Toeplitz-plus-finite-rank Jacobi operators like $J$ and $D$ occurring in Theorem 4.8 and Corollary 4.10, the operators $C$, $C^{T}$, $C^{\mu}$ and $(C^{\mu})^{T}$ all map $\ell_{\mathcal{F}}$ to $\ell_{\mathcal{F}}$. Consequently, combinations such as $CC^{T}$, $C^{\mu}C^{T}$ are all well defined operators from $\ell_{\mathcal{F}}$ to $\ell_{\mathcal{F}}$.

4.4 Properties of the resolvent

When the Jacobi operator $J$ is Toeplitz-plus-finite rank, as a consequence of the structure of the connection coefficients matrix proved in subsection 4.3, the principal resolvent $G$ (see Definition 2.1) and spectral measure (see Theorem 2.2) are also highly structured. As usual these proofs are stated for a finite-rank perturbation of the free Jacobi operator $\Delta$, but apply to general Toeplitz-plus-finite rank Jacobi operators by applying appropriate scaling and shifting.

Theorem 4.12**.**

Let $J$ be a Jacobi operator such that there exists an $n$ such that

[TABLE]

i.e. it is equal to the free Jacobi operator $\Delta$ outside the $n\times n$ principal submatrix. Then the principal resolvent for $J$ is

[TABLE]

where

[TABLE]

$P_{k}$* are the orthonormal polynomials for $J$, $P_{k}^{\mu}$ are the first associated polynomials for $J$ as in Definition 2.5, and $U_{k}$ are the Chebyshev polynomials of the second kind.*

Remark 4.13*.*

$p_{C}^{\mu}$ is the $\mu$-derivative of $p_{C}$ as in Definition 2.7.

Proof.

Using Theorem 2.2 and Proposition 3.6,

[TABLE]

Now, since $p_{C}$ is a polynomial we can split this into

[TABLE]

The first term is equal to $p_{C}(\lambda)G(\lambda)$, and the second term is equal to $p_{C}^{\mu}(\lambda)$ by Lemma 2.6 and Remark 4.13. The equation can now be immediately rearranged to obtain (4.9).

To see the equality in equation (4.10), note that by the definition of the connection coefficient matrix $C$,

[TABLE]

Equation (4.11) follows by the same algebra. ∎

Theorem 4.14**.**

Let $J$ be a Jacobi operator such that there exists an $n$ such that

[TABLE]

i.e. it is equal to the free Jacobi operator $\Delta$ outside the $n\times n$ principal submatrix. Then the spectral measure for $J$ is

[TABLE]

where $\lambda_{1},\ldots,\lambda_{r}$ are the roots of $p_{C}$ in $\mathbb{R}\setminus\{1,-1\}$ such that

[TABLE]

There are no roots of $p_{C}$ inside $(-1,1)$, but there may be simple roots at $\pm 1$.

Remark 4.15*.*

We will see in Theorem 4.22 that the number of roots of $p_{C}$ for which $w_{k}\neq 0$ is at most $n$ (i.e. $r\leq n$). Hence, while the degree of $p_{C}$ is at most $2n-1$, at least $n-1$ are cancelled out by factors in the numerator.

Proof.

Let $G$ and $\mu$ be the principal resolvent and spectral measure of $J$ respectively. By Theorem 4.12,

[TABLE]

Letting $\lambda_{1},\ldots,\lambda_{2n-1}$ be the roots of $p_{C}$ in the complex plane, define the set

[TABLE]

By inspection of the above formula for $G$, and because resolvents of selfadjoint operators are analytic off the real line, we have that $G$ is continuous outside of $S$. Therefore, for any $s\in\mathbb{R}$ such that $\mathrm{dist}(s,S)>0$, we have

[TABLE]

Hence by Theorem 2.2 part (ii), for any interval $(s_{1},s_{2})$ such that $\mathrm{dist}(S,(s_{1},s_{2}))>0$, we have $\mu((s_{1},s_{2}))+\frac{1}{2}\mu(\{s_{1}\})+\frac{1}{2}\mu(\{s_{2}\})=0$. Therefore the essential support of $\mu$ is contained within $S$.

We are interested in the real roots of $p_{C}$. Let us consider the potential for roots of $p_{C}$ in the interval $[-1,1]$. By Proposition 3.6, $\mathrm{d}\mu_{\Delta}(s)=p_{C}(s)\mathrm{d}\mu(s)$ for all $s\in\mathbb{R}$. For any $s\in[-1,1]$ such that $p_{C}(s)\neq 0$, it follows that $\mathrm{d}\mu(s)=\frac{2}{\pi}\frac{\sqrt{1-s^{2}}}{p_{C}(s)}\mathrm{d}s$. From this we have

[TABLE]

This integral is only finite, so $p_{C}$ has no roots in $(-1,1)$, but may have simple roots at $\pm 1$. One example where we have simple roots at $\pm 1$ is seen in Example 4.18 with $\beta=\sqrt{2}$.

Since $S$ is a disjoint union of $[-1,1]$ and a finite set $S^{\prime}$ we can write

[TABLE]

By Theorem 2.2 part (iii), $\mu(\{s\})=\lim_{\epsilon\searrow 0}\frac{\epsilon}{i}G(s+i\epsilon)\text{ for all }s\in\mathbb{R}$. This gives the desired formula for $w_{k}$. ∎

Remark 4.16*.*

Theorem 4.14 gives an explicit formula for the spectral measure of $J$, when $J$ is Toeplitz-plus-finite-rank Jacobi operator. The entries of $C$ can be computed in $\mathcal{O}(n^{2})$ operations (for an $n\times n$ perturbation of Toeplitz). Hence, the absolutely continuous part of the measure can be computed exactly in finite time. It would appear at first that we may compute the locations of the point spectrum by computing the roots of $p_{C}$, but as stated in Remark 4.15 we find that not all real roots of $p_{C}$ have $w_{k}\neq 0$. Hence we rely on cancellation between the numerator and denominator in the formula for $G(\lambda)$, which is a dangerous game, because if roots of polynomials are only known approximately then it is impossible to distinguish between cancellation and the case where a pole and a root are merely extremely close. Subsection 4.5 remedies this situation.

Example 4.17** (Basic perturbation 1 revisited).**

The polynomial $p_{C}$ in Theorem 4.12 is

[TABLE]

and the $\mu$-derivative is $p_{C}^{\mu}(\lambda)=-2\alpha$. Theorem 4.12 gives

[TABLE]

Consider the case $|\alpha|\leq 1$. Then a brief calculation reveals $G_{\Delta}(\frac{1}{2}(\alpha+\alpha^{-1}))=-2\alpha$. Hence the root $\lambda=\frac{1}{2}(\alpha+\alpha^{-1})$ of the denominator is always cancelled out. Hence $G$ has no poles, and so $J$ has no eigenvalues.

In the case where $|\alpha|>1$, we have a different situation. Here $G_{\Delta}(\frac{1}{2}(\alpha+\alpha^{-1}))=-2\alpha^{-1}$. Therefore the root $\lambda=\frac{1}{2}(\alpha+\alpha^{-1})$ of the denominator is never cancelled out. Hence there is always a pole of $G$ at $\lambda=\frac{1}{2}(\alpha+\alpha^{-1})$, and therefore also an eigenvalue of $J$ there.

Notice a heavy reliance on cancellations in the numerator and denominator for the existence of eigenvalues. The approach in subsection 4.5 avoids this.

Example 4.18** (Basic perturbation 2 revisited).**

The polynomial $p_{C}$ in Theorem 4.12 is

[TABLE]

This simplifies to $p_{C}(\lambda)=4(1-\beta^{-2})\left(\frac{\beta^{4}}{4(\beta^{2}-1)}-\lambda^{2}\right)$. Using Definition 2.5, the $\mu$-derivative is $p_{C}^{\mu}(\lambda)=c_{0,2}P_{2}^{\mu}(\lambda)=4\beta^{-1}\lambda$. Theorem 4.12 gives

[TABLE]

Clearly the only points at which $G$ may have a pole is $\lambda=\pm\frac{\beta^{2}}{2\sqrt{\beta^{2}-1}}$. However, it is difficult to see whether there would be cancellation on the numerator. In the previous discussion on this example we noted that there would not be any poles when $|\beta|<\sqrt{2}$, which means that the numerator must be zero at these points, but it is far from clear here. The techniques we develop in the sequel illuminate this issue, especially for examples which are much more complicated than the two trivial ones given so far.

4.5 The Joukowski transformation

The following two lemmata and two theorems prove that the issue of cancellation and the number of discrete spectra in Theorem 4.12 and Theorem 4.14 can be solved by making the change of variables

[TABLE]

This map is known as the Joukowski map. It is an analytic bijection from $\mathbb{D}=\{z\in\mathbb{C}:|z|<1\}$ to $\mathbb{C}\setminus[-1,1]$, sending the unit circle to two copies of the interval $[-1,1]$.

The Joukowski map has special relevance for the principal resolvent of $\Delta$. A brief calculation reveals that for $z\in\mathbb{D}$,

[TABLE]

Further, we will see that the polynomials $p_{C}(\lambda)$ and $p_{C}^{\mu}(\lambda)$ occurring in our formula for $G$ can be expressed neatly as polynomials in $z$ and $z^{-1}$. This is a consequence of a special property of the Chebyshev polynomials of the second kind, that for any $k\in\mathbb{Z}$ and $z\in\mathbb{D}$

[TABLE]

These convenient facts allow us to remove any square roots involved in the formulae in Theorem 4.12.

Lemma 4.19**.**

Let $p_{C}(\lambda)=\sum_{k=0}^{2n-1}\langle e_{0},CC^{T}e_{k}\rangle U_{k}(\lambda)$ as in Theorem 4.12 and let $c$ be the symbol of $C_{\rm Toe}$, the Toeplitz part of $C$ as guaranteed by Theorem 4.8. Then

[TABLE]

where $\lambda(z)=\frac{1}{2}(z+z^{-1})$.

Proof.

The key quantity to observe for this proof is

[TABLE]

for $z\in\mathbb{D}$ as $m\to\infty$. We will show it is equal to both sides of equation (4.15). Consider the polynomial $U_{m}\cdot p_{C}$. The $j$th coefficient in an expansion in the basis $P_{0},P_{1},P_{2},\ldots$ is given by

[TABLE]

because $p_{C}=\frac{\mathrm{d}\mu_{\Delta}}{\mathrm{d}\mu}$ by Proposition 3.6. Hence

[TABLE]

for all $m\in\mathbb{N}$ and all $z\in\mathbb{D}$.

Now we show that (4.16) converges to $c(z)c(z^{-1})$ as $m\to\infty$. By the definition of the connection coefficients, $P_{m+k}=\sum_{j=0}^{2n-1}c_{m+k-j,m+k}U_{m+k-j}$. Therefore,

[TABLE]

Now, by Theorem 4.8, $C=C_{\rm Toe}+C_{\rm fin}$, where $C_{\rm fin}$ is zero outside the $(n-1)\times(2n-2)$ principal submatrix. Hence for $m$ sufficiently large we have $c_{m,m+k}=t_{k}$ for a sequence $(t_{k})_{k\in\mathbb{Z}}$ such that $t_{k}=0$ for $k\notin\{0,1,\ldots,2n-1\}$. Hence we have for $m$ sufficiently large,

[TABLE]

By equation (4.14), this tends to $\sum_{k=0}^{2n-1}\sum_{j=0}^{2n-1}t_{k}t_{j}z^{j-k}$ as $m\to\infty$. This is equal to $c(z)c(z^{-1})$, as required to complete the proof. ∎

Lemma 4.20**.**

Let $p^{\mu}_{C}(\lambda)=\sum_{k=0}^{2n-1}\langle e_{k},C^{\mu}C^{T}e_{0}\rangle U_{k}(\lambda)$ as in Theorem 4.12 and let $c_{\mu}$ be the symbol of $C_{\rm Toe}^{\mu}$, the Toeplitz part of $C^{\mu}$ as guaranteed by Corollary 4.10. Then

[TABLE]

where $\lambda(z)=\frac{1}{2}(z+z^{-1})$ and $z\in\mathbb{D}$.

Proof.

The key quantity to observe for this proof is

[TABLE]

for $z\in\mathbb{D}$, as $m\to\infty$. We will compute two equivalent expressions for this quantity to derive equation (4.17). In the proof of Lemma 4.19, it was shown that $U_{m}(\lambda)p_{C}(\lambda)=\sum_{k=0}^{2n-1}c_{m,m+k}P_{k}(\lambda)$. We take the $\mu$-derivative (see Definition 2.7) of both sides as follows.

[TABLE]

because by Proposition 3.6, $p_{C}=\frac{\mathrm{d}\mu_{\Delta}}{\mathrm{d}\mu}$. Using the formula from equation (4.3), $U_{m}^{\mu_{\Delta}}=2U_{m-1}$, we find that the $\mu$-derivative of $U_{m}(s)p_{C}(s)$ is equal to $U_{m}(\lambda)p^{\mu}_{C}(\lambda)+2U_{m-1}(\lambda)$. Taking the $\mu$-derivative from the other side gives $\sum_{k=0}^{2n-1}c_{m,m+k}P^{\mu}_{k}(\lambda)$. Taking the limit as $m\to\infty$ and using equation (4.14), we have our first limit for the quantity in equation (4.18):

[TABLE]

Now we show that (4.18) converges to $c_{\mu}(z)c(z^{-1})$ as $m\to\infty$. By the definition of the connection coefficients matrix $C^{\mu}$, $P^{\mu}_{m+k}=\sum_{j=1}^{2n-2}c^{\mu}_{m+k-j,m+k}U_{m+k-j}$. Therefore,

[TABLE]

By Corollary 4.10, $C^{\mu}=C^{\mu}_{\rm Toe}+C^{\mu}_{\rm fin}$, where $C^{\mu}_{\rm fin}$ is zero outside the principal $(n-2)\times(2n-1)$ submatrix. Hence for $m$ sufficiently large we have $c^{\mu}_{m,m+k}=t^{\mu}_{k}$ for a sequence $(t^{\mu}_{k})_{k\in\mathbb{Z}}$ such that $t^{\mu}_{k}=0$ for $k\notin\{1,\ldots,2n-2\}$. Hence we have for sufficiently large $m$,

[TABLE]

By equation (4.14), this tends to $\sum_{k=0}^{2n-1}\sum_{j=1}^{2n-2}t_{k}t^{\mu}_{j}z^{j-k}$ as $m\to\infty$. This is equal to $c_{\mu}(z)c(z^{-1})$. Equating this with the other equality for equation (4.18) gives $p^{\mu}(\lambda(z))=c(z^{-1})c_{\mu}(z)-2z$ as required. ∎

The following theorem describes Theorem 4.12 under the change of variables induced by the Joukowski map. The remarkable thing is that the resolvent is expressible as a rational function inside the unit disc.

Theorem 4.21**.**

Let $J$ be a Jacobi operator such that there exists an $n$ such that

[TABLE]

i.e. it is equal to the free Jacobi operator $\Delta$ outside the $n\times n$ principal submatrix. By Theorem 4.8 the connection coefficient matrix can be decomposed into $C=C_{\rm Toe}+C_{\rm fin}$. By Corollary 4.10, we similarly have $C^{\mu}=C^{\mu}_{\rm Toe}+C^{\mu}_{\rm fin}$. If $c$ and $c_{\mu}$ are the Toeplitz symbols of $C_{\rm Toe}$ and $C^{\mu}_{\rm Toe}$ respectively, then for $\lambda(z)=\frac{1}{2}(z+z^{-1})$ with $z\in\mathbb{D}$, the principal resolvent $G$ is given by the rational function

[TABLE]

Proof.

Combining Theorem 4.12, equation (4.13) and Lemmata 4.19 and 4.20, we have

[TABLE]

This completes the proof. ∎

The following theorem gives a better description of the weights $w_{k}$ in Theorem 4.14, utilising the Joukowski map and the Toeplitz symbol $c$.

Theorem 4.22**.**

Let $J$ be a Jacobi operator such that there exists an $n$ such that

[TABLE]

i.e. it is equal to the free Jacobi operator $\Delta$ outside the $n\times n$ principal submatrix. By Theorem 4.8 the connection coefficient matrix can be written $C=C_{\rm Toe}+C_{\rm fin}$. If $c$ is the Toeplitz symbol of $C_{\rm Toe}$, then the spectral measure of $J$ is

[TABLE]

Here $z_{k}$ are the roots of $c$ that lie in the open unit disk, which are all real and simple. The only roots of $c$ on the unit circle are $\pm 1$, which can also only be simple. Further, $r\leq n$.

Proof.

By Theorem 4.14,

[TABLE]

where $r\leq n$. Hence we just need to prove something more specific about the roots of $c$, $\lambda_{1},\ldots,\lambda_{r}$, and $w_{1},\ldots,w_{r}$.

By Theorem 4.21, $G(\lambda(z))=-c_{\mu}(z)/c(z)$ for $z\in\mathbb{D}$. By Lemma 4.20, $c(z^{-1})c_{\mu}(z)-2z=p^{\mu}(\lambda(z))=p^{\mu}(\lambda(z^{-1}))=c(z)c_{\mu}(z^{-1})-2z^{-1}$, so

[TABLE]

Therefore $c(z)$ and $c_{\mu}(z)$ cannot simultaneously be zero unless $z=z^{-1}$, which only happens at $z=\pm 1$. By the same reasoning, $c(z)$ and $c(z^{-1})$ also cannot be simultaneously zero unless $z=\pm 1$. Since the Joukowski map $\lambda$ is a bijection from $\mathbb{D}$ to $\mathbb{C}\setminus[-1,1]$, this shows that the (simple and real) poles of $G$ in $\mathbb{C}\setminus[-1,1]$ are precisely $\lambda(z_{1}),\ldots,\lambda(z_{r})$, where $z_{1},\ldots,z_{r}$ are the (necessarily simple and real) roots of $c$ in $\mathbb{D}$.

What are the values of the weights of the Dirac deltas, $w_{1},\ldots,w_{r}$? By Theorem 4.14,

[TABLE]

By equation (4.20), since $c(z_{k})=0$, we have $c_{\mu}(z_{k})=2(z_{k}-z_{k}^{-1})/c(z_{k}^{-1})$. This gives

[TABLE]

Note that if $c(z)=0$ then $c(\overline{z})=0$ because $c$ has real coefficients. If $c$ has a root $z_{0}$ on the unit circle, then $c(z_{0})=c(z_{0}^{-1})=0$ because $\overline{z_{0}}=z_{0}^{-1}$, which earlier in the proof we showed only occurs if $z_{0}=\pm 1$. Hence $c$ does not have roots on the unit circle except possibly $\pm 1$. ∎

Example 4.23** (Basic perturbation 1 re-revisited).**

Considering the connection coefficient matrix in equation (4.4), we see that the Toeplitz symbol $c$ is $c(z)=1-\alpha z$. By Theorem 4.22 the roots of $c$ in the unit disc correspond to eigenvalues of $J_{\alpha}$. As is consistent with our previous considerations, $c$ has a root in the unit disc if and only if $|\alpha|>1$, and those eigenvalues are $\lambda(\alpha^{-1})=\frac{1}{2}(\alpha+\alpha^{-1})$. See Appendix A for figures depicting the spectral measure and the resolvent.

Example 4.24** (Basic perturbation 2 re-revisited).**

Considering the connection coefficient matrix in equation (4.5), we see that the Toeplitz symbol $c$ is $c(z)=\beta^{-1}+(\beta^{-1}-\beta)z^{2}$. By Theorem 4.22 the roots of $c$ in the unit disc correspond to eigenvalues of $J_{\beta}$. The roots of $c$ are $\pm\frac{1}{\sqrt{\beta^{2}-1}}$. If $\beta\in\left(0,\sqrt{2}\right]\setminus\{1\}$ then $\left|\pm\frac{1}{\sqrt{\beta^{2}-1}}\right|\geq 1$ so there are no roots of $c$ in the unit disc, as is consistent with the previous observations. What was difficult to see before is, if $\beta>\sqrt{2}$ then $\left|\pm\frac{1}{\sqrt{\beta^{2}-1}}\right|<1$, so there is a root of $c$ inside $\mathbb{D}$, and it corresponds to an eigenvalue,

[TABLE]

See Appendix A for figures depicting the spectral measure and the resolvent.

5 Toeplitz-plus-trace-class Jacobi operators

In this section we extend the results of the previous section to the case where the Jacobi operator is Toeplitz-plus-trace-class. This cannot be done as a direct extension of the work in the previous section as the formulae obtained depended on the fact that some of the functions involved were merely polynomials in order to have a function defined for all $\lambda$ in an a priori known region of the complex plane. We admit that it may be possible to perform the analysis directly, but state that it is not straightforward. We are interested in feasible (finite) computation so are content to deal directly with the Toeplitz-plus-finite-rank case and perform a limiting process. The crucial question for computation is, can we approximate the spectral measure of a Toeplitz-plus-trace-class Jacobi operator whilst reading only finitely many entries of the matrix?

Here we make clear the definition of a Toeplitz-plus-trace-class Jacobi operator.

Definition 5.1**.**

An operator $K:\ell^{2}\to\ell^{2}$ is said to be trace class if $\sum_{k=0}^{\infty}e_{k}^{T}(K^{T}K)^{1/2}e_{k}<\infty$. Hence we say that a Jacobi operator $J$ such that $\alpha_{k}\to 0$, $\beta_{k}\to\frac{1}{2}$ as $k\to\infty$ is Toeplitz-plus-trace-class if

[TABLE]

5.1 Jacobi operators for Jacobi polynomials

The most well known class of orthogonal polynomials is the Jacobi polynomials, whose measure of orthogonality is

[TABLE]

where $\alpha$,$\beta>-1$ and $B$ is Euler’s Beta function. The Jacobi operator for the normalised Jacobi polynomials with respect to this probability measure, and hence the three-term recurrence coefficients, are given by [36],

[TABLE]

Note that $|\alpha_{k}|=\mathcal{O}(k^{-2})$ and

[TABLE]

Hence the Jacobi operators for the Jacobi polynomials are Toeplitz-plus-trace-class for all $\alpha,\beta>-1$.

The Chebyshev polynomials $T_{k}$ and $U_{k}$ discussed in the previous section are specific cases of Jacobi polynomials, with $\alpha,\beta=-\frac{1}{2},-\frac{1}{2}$ for $T_{k}$ and $\alpha,\beta=\frac{1}{2},\frac{1}{2}$ for $U_{k}$.

In Appendix A numerical computations of the spectral measures and resolvents of these Jacobi operators are presented.

5.2 Toeplitz-plus-finite-rank approximations

We propose to use the techniques from Section 4. Therefore for a Jacobi operator $J$, we can define the Toeplitz-plus-finite-rank approximations $J^{[m]}$, where

[TABLE]

Each Jacobi operator $J^{[m]}$ has a spectral measure $\mu^{[m]}$ which can be computed using Theorem 4.22. The main question for this section is: how do the computable measures $\mu^{[m]}$ approximate the spectral measure $\mu$ of $J$?

Proposition 5.2**.**

Let $J$ a Jacobi operator (bounded, but with no assumed structure imposed) and let $\mu$ be its spectral measure. Then the measures $\mu^{[1]},\mu^{[2]},\ldots$ which are the spectral measures of $J^{[1]},J^{[2]},\ldots$ converge to $\mu$ in a weak sense. Precisely,

[TABLE]

for all $f\in C_{b}(\mathbb{R})$.

Proof.

Each spectral measure $\mu^{[m]}$ and $\mu$ are supported on the spectra of $J^{[m]}$ and $J$, each of which are contained within $[-\|J^{[m]}\|_{2},\|J^{[m]}\|_{2}]$ and $[-\|J\|_{2},\|J\|_{2}]$. Since $\|J^{[m]}\|_{2}$ and $\|J\|_{2}$ are less than

[TABLE]

we have that all the spectral measures involved are supported within the interval $[-M,M]$. Hence we can consider integrating functions $f\in C([-M,M])$ without ambiguity.

By Weierstrass’ Theorem, polynomials are dense in $C([-M,M])$, so we only need to consider polynomials as test functions, and by linearity we only need to consider the orthogonal polynomials for $J$. The first polynomial $P_{0}$ has immediate convergence, since the measures are all probability measures. Now consider $P_{k}$ for some $k>0$, which satisfies $\int P_{k}(s)\,\mathrm{d}\mu(s)=0$. For $m>k$, $P_{k}$ is also the $k$th orthogonal polynomial for $J^{[m]}$, hence $\int P_{k}(s)\,\mathrm{d}\mu^{[m]}(s)=0$. This completes the proof. ∎

5.3 Asymptotics of the connection coefficients

Here we formulate a lower triangular block operator equation $\mathcal{L}\underline{c}=e_{0}^{0}$ , where $e_{0}^{0}=(e_{0},0,0,\ldots)^{\top}$, satisfied by the entries of the connection coefficient matrices encoded into a vector $\underline{c}$. For Toeplitz-plus-trace-class Jacobi operators we give appropriate Banach spaces upon which the operator $\mathcal{L}$ is bounded and invertible, enabling precise results about the asymptotics of the connection coefficients to be derived.

Lemma 5.3**.**

Let $J$ and $D$ be Jacobi operators with entries $\{\alpha_{k},\beta_{k}\}_{k=0}^{\infty}$ and $\{\gamma_{k},\delta_{k}\}_{k=0}^{\infty}$ respectively. If we decompose the upper triangular part of $C_{J\to D}$ into a sequence of sequences, stacking each diagonal on top of each other, we get the following block linear system,

[TABLE]

where for each $i$,

[TABLE]

For $B_{-1}$ to make sense we define $\beta_{-1}=1/2$.

Proof.

This is simply the 5-point discrete system in Lemma 3.2 rewritten. ∎

We write the infinite-dimensional-block-infinite-dimensional system (5.2) in the form,

[TABLE]

For general Jacobi operators $J$ and $D$, the operators $A_{i}$ and $B_{i}$ are well defined linear operators from $\ell_{\mathcal{F}}^{\star}$ to $\ell_{\mathcal{F}}^{\star}$. The block operator $\mathcal{L}$ is whence considered as a linear operator from the space of sequences of real sequences, $\ell_{\mathcal{F}}^{\star}(\ell_{\mathcal{F}}^{\star})$ to itself. We will use this kind of notation for other spaces as follows.

Definition 5.4** (Vector-valued sequences).**

If $\ell_{X}$ is a vector space of scalar-valued sequences, and $Y$ is another vector space then we let $\ell_{X}(Y)$ denote the vector space of sequences of elements of $Y$. In many cases in which $\ell_{X}$ and $Y$ are both normed spaces, then $\ell_{X}(Y)$ naturally defines a normed space in which the norm is derived from that of $\ell_{X}$ by replacing all instances of absolute value with the norm on $Y$. For example, $\ell^{p}(\ell^{\infty})$ is a normed space with norm $\|(a_{k})_{k=0}^{\infty}\|_{\ell^{p}(\ell^{\infty})}=\left(\sum_{k=0}^{\infty}\|a_{k}\|_{\infty}^{p}\right)^{\frac{1}{p}}$.

The following two spaces are relevant for the Toeplitz-plus-trace-class Jacobi operators.

Definition 5.5** (Sequences of bounded variation).**

Following [21, Ch. IV.2.3], denote by $bv$ the Banach space of all sequences with bounded variation, that is sequences such that the norm

[TABLE]

is finite.

The following result is immediate from the definition of the norm on $bv$.

Lemma 5.6**.**

*There is a continuous embedding of $bv$ into the Banach space of convergent sequences (endowed with the supremum norm) i.e. for all $(a_{k})_{k=0}^{\infty}\in bv$, $\lim_{k\to\infty}a_{k}$ exists, and $\sup_{k}|a_{k}|\leq\|(a_{k})_{k=0}^{\infty}\|_{bv}$. Furthermore, $\lim_{k\to\infty}|a_{k}|\leq\|a\|_{bv}$. *

Definition 5.7** (Geometrically weighted $\ell^{1}$).**

For any $R>0$, e define the Banach space $\ell^{1}_{R}$ to be the space of sequences such that the norm

[TABLE]

is finite.

Proposition 5.8**.**

The operator norm on $\ell^{1}_{R}$ is equal to

[TABLE]

The following Lemma and its Corollary show that it is natural to think of $\underline{c}$ as lying in the space $\ell^{1}_{R}(bv)$.

Lemma 5.9**.**

Let $J=\Delta+K$ be a Jacobi operator where $K$ is trace class and let $D=\Delta$. Then for any $R\in(0,1)$ the operator $\mathcal{L}$ in equation (5.3) is bounded and invertible as an operator from $\ell_{R}^{1}(bv)$ to $\ell^{1}_{R}(\ell^{1})$. Furthermore, if $\mathcal{L}^{[m]}$ is the operator in equation (5.3) generated by the Toeplitz-plus-finite-rank truncation $J^{[m]}$, then

[TABLE]

Proof.

We can write $\mathcal{L}$ in equation (5.3) in the form $\mathcal{L}=\mathcal{T}+\mathcal{K}$ where

[TABLE]

and

[TABLE]

This decomposition will allow us to prove that $\mathcal{L}$ is bounded and invertible as follows. We will show that as operators from $\ell^{1}_{R}(bv)$ to $\ell^{1}_{R}(\ell^{1})$, $\mathcal{T}$ is bounded and invertible, and $\mathcal{K}$ is compact. This implies that $\mathcal{L}$ is a Fredholm operator with index [math]. Therefore, by the Fredholm Alternative Theorem, $\mathcal{L}$ is invertible if and only if it is injective. It is indeed injective, because it is block lower triangular with invertible diagonal blocks, so forward substitution on the system $\mathcal{L}\underline{v}=\underline{0}$ implies that each entry of $\underline{v}$ must be zero.

First let us prove that $\mathcal{T}$ is bounded and invertible. It is elementary that $T$ is an isometric isomorphism from $bv$ to $\ell^{1}$ and $T^{T}$ is bounded with norm at most 1. Hence using Proposition 5.8 we have

[TABLE]

Because each operator is lower triangular, the left and right inverse of $\mathcal{T}:\ell_{\mathcal{F}}(\ell_{\mathcal{F}})\rightarrow\ell_{\mathcal{F}}(\ell_{\mathcal{F}})$ is

[TABLE]

This matrix is block-lower triangular and block-Toeplitz with first column having $2i$th block of the form $T^{-1}(-T^{T}T^{-1})^{i}$ and $(2i+1)$th block zero. We must check that $\mathcal{T}^{-1}$ is bounded in the norms on $\ell_{R}^{1}(\ell^{1})$ to $\ell_{R}^{1}(bv)$ so that it may be extended to $\ell_{R}^{1}(\ell^{1})$ from the dense subspace $\ell_{\mathcal{F}}$. Again using Proposition 5.8 we have

[TABLE]

Now let us prove that $\mathcal{K}:\ell_{R}^{1}(bv)\to\ell_{R}^{1}(\ell^{1})$ is compact. Consider the finite rank operator $\mathcal{K}^{[m]}$, where all elements are the same as in $\mathcal{K}$, except that all occurrences of $\alpha_{i}$ and $2\beta_{i}-1$ are replaced by [math] for $i\geq m$. Using Proposition 5.8 we have

[TABLE]

By the continuous embedding in Lemma 5.6, $\|\cdot\|_{bv\to\ell^{1}}\leq\|\cdot\|_{\ell^{\infty}\to\ell^{1}}$. Hence

[TABLE]

This convergence is due to the fact that $J-\Delta$ is trace class. Since $\mathcal{K}$ is a norm limit of finite rank operators it is compact. This completes the proof that $\mathcal{L}$ is bounded and invertible.

Now consider the operator $\mathcal{L}^{[m]}$ defined in the statement of the Lemma, which is equal to $\mathcal{T}+\mathcal{K}^{[m]}$ (where $\mathcal{K}^{[m]}$ is precisely that which was considered whilst proving $\mathcal{K}$ is compact). Hence,

[TABLE]

This completes the proof. ∎

Corollary 5.10**.**

Let $J=\Delta+K$ be a Jacobi operator where $K$ is trace class and let $\underline{c}\in\ell_{\mathcal{F}}^{\star}(\ell_{\mathcal{F}}^{\star})$ be the vector of diagonals of $C_{J\to\Delta}$ as in equation (5.3). Then $\underline{c}\in\ell^{1}_{R}(bv)$. If $J$ has Toeplitz-plus-finite-rank approximations $J^{[m]}$ and $\underline{c}^{[m]}$ denotes the vector of diagonals of $C^{[m]}$, then

[TABLE]

Proof.

By equation (5.3)

[TABLE]

Since $\|e_{0}^{0}\|_{\ell_{R}^{1}(\ell^{1})}=1$, the proof is completed if we show $\|\mathcal{L}^{-1}-(\mathcal{L}^{[m]})^{-1}\|_{\ell_{R}^{1}(\ell^{1})\to\ell_{R}^{1}(bv)}\to 0$ as $m\to\infty$.

Suppose that $m$ is sufficiently large so that $\|\mathcal{L}-\mathcal{L}^{[m]}\|<\|\mathcal{L}^{-1}\|^{-1}$ (guaranteed by Lemma 5.9). Note that $\mathcal{L}^{-1}$ is bounded by the Inverse Mapping Theorem and Lemma 5.9. Then by a well-known result (see for example, [2]),

[TABLE]

This tends to zero as $m\to\infty$, by Lemma 5.9. ∎

Theorem 5.11**.**

Let $J=\Delta+K$ be a Jacobi operator where $K$ is trace class. Then $C=C_{J\to\Delta}$ can be decomposed into

[TABLE]

where $C_{\rm Toe}$ is upper triangular, Toeplitz and bounded as an operator from $\ell^{1}_{R}$ to $\ell^{1}_{R}$, and $C_{\rm com}$ is compact as an operator from $\ell^{1}_{R}$ to $\ell^{1}_{R}$, for all $R>1$. Also, if $J$ has Toeplitz-plus-finite-rank approximations $J^{[m]}$ with connection coefficient matrices $C^{[m]}=C_{\rm Toe}^{[m]}+C_{\rm com}^{[m]}$, then

[TABLE]

in the operator norm topology over $\ell^{1}_{R}$ for all $R>1$.

Proof.

By Lemma 5.9, for each $k$ the sequence $(c_{0,0+k},c_{1,1+k},c_{2,2+k},\ldots)$ is an element of $bv$. By Lemma 5.6 each is therefore a convergent sequence, whose limits we call $t_{k}$. Hence we can define an upper triangular Toeplitz matrix $C_{\rm Toe}$ whose $(i,j)$th element is $t_{j-i}$, and define $C_{\rm com}=C-C_{\rm Toe}$.

The Toeplitz matrix $C_{\rm Toe}$ is bounded from $\ell^{1}_{R}$ to $\ell^{1}_{R}$ for all $R>1$ by the following calculation.

[TABLE]

By Lemma 5.9 this quantity is finite (since $R^{-1}\in(0,1)$).

Now we show convergence results. The compactness of $C_{\rm com}$ will follow at the end. For all $R>1$,

[TABLE]

For the third line of the above sequence of equations, note that for fixed $k$, $c_{0,k}-c^{[m]}_{0,k},c_{1,1+k}-c^{[m]}_{1,1+k},c_{2,2+k}-c^{[m]}_{2,2+k},\ldots$ is a $bv$ sequence, and refer to Lemma 5.6.

[TABLE]

For the third line of the above sequence, note that $t_{k}-t_{k}^{[m]}$ is the limit of the $bv$ sequence $c_{*,*+k}-c^{[m]}_{*,*+k}$, and refer to Lemma 5.6.

[TABLE]

Using Corollary 5.10, that $\|\underline{c}-\underline{c}^{[m]}\|_{\ell^{1}_{R}(bv)}\to 0$ as $m\to\infty$, we have the convergence results.

By Theorem 4.8, $C_{\rm com}^{[m]}$ has finite rank. Therefore, since $C_{\rm com}=\lim_{m\to\infty}C_{\rm com}^{[m]}$ in the operator norm topology over $\ell^{1}_{R^{-1}}$, we have that $C_{\rm com}$ is compact in that topology. ∎

Corollary 5.12**.**

Let $C^{\mu}$ be as defined in Definition 3.11 for $C$ as in Theorem 5.11. Then $C^{\mu}$ can be decomposed into $C^{\mu}=C^{\mu}_{\rm Toe}+C^{\mu}_{\rm com}$ where $C^{\mu}_{\rm Toe}$ is upper triangular, Toeplitz and bounded as an operator from $\ell^{1}_{R}$ to $\ell^{1}_{R}$, and $C^{\mu}_{\rm com}$ is compact as an operator from $\ell^{1}_{R}$ to $\ell^{1}_{R}$, for all $R>1$. Furthermore, if $J$ has Toeplitz-plus-finite-rank approximations $J^{[m]}$ with connection coefficient matrices $(C^{\mu})^{[m]}=(C^{\mu}_{\rm Toe})^{[m]}+(C^{\mu}_{\rm com})^{[m]}$, then

[TABLE]

in the operator norm topology over $\ell^{1}_{R}$.

Proof.

This follows from Theorem 5.11 applied to $J^{\mu}$ as defined in Lemma 3.13. ∎

Theorem 5.13**.**

Let $J$ be a Jacobi operator such that $J=\Delta+K$ where $K$ is trace class. The Toeplitz symbols $c$ and $c_{\mu}$ of the Toeplitz parts of $C_{J\to\Delta}$ and $C^{\mu}_{J\to\Delta}$ are both analytic in the unit disc. Furthermore, if $J$ has Toeplitz-plus-finite-rank approximations $J^{[m]}$ with Toeplitz symbols $c^{[m]}$ and $c_{\mu}^{[m]}$, then $c^{[m]}\to c$ and $c^{[m]}_{\mu}\to c_{\mu}$ as $m\to\infty$ uniformly on compact subsets of $\mathbb{D}$.

Proof.

Let $R>1$, and let $0\leq|z|\leq R^{-1}<1$. Then by Lemma 5.6 we have

[TABLE]

where $\underline{c}$ is as defined in equation (5.3). By Lemma 5.9 this quantity is finite. Since $R$ is arbitrary, the radius of convergence of the series is 1. The same is true for $c_{\mu}$ by Lemma 3.13.

Now we prove that the Toeplitz symbols corresponding to the Toeplitz-plus-finite-rank approximations converge.

[TABLE]

To go between the first and second lines, note that for each $k$, $c_{*,*+k}-c^{[m]}_{*,*+k}$ is a $bv$ sequence whose limit is $t_{k}-t^{[m]}_{k}$ and refer to Lemma 5.6. Now, $\|\underline{c}-\underline{c}^{[m]}\|_{\ell_{R^{-1}}^{1}(bv)}\to 0$ as $m\to\infty$ by Corollary 5.10. The same is true for $\sup_{|z|\leq R^{-1}}|c_{\mu}(z)-c_{\mu}^{[m]}(z)|$ by Lemma 3.13. ∎

Theorem 5.14** (See [31]).**

Let $A$ and $B$ be bounded self-adjoint operators on $\ell^{2}$. Then

[TABLE]

Theorem 5.15**.**

Let $J=\Delta+K$ be a Toeplitz-plus-trace-class Jacobi operator, and let $c$ and $c_{\mu}$ be the analytic functions as defined in Theorem 5.11 and Corollary 5.12. Then for $\lambda(z)=\frac{1}{2}(z+z^{-1})$ with $z\in\mathbb{D}$ such that $\lambda(z)\notin\sigma(J)$, the principal resolvent $G$ is given by the meromorphic function

[TABLE]

Therefore, all eigenvalues of $J$ are of the form $\lambda(z_{k})$, where $z_{k}$ is a root of $c$ in $\mathbb{D}$.

Proof.

Let $z\in\mathbb{D}$ such that $\lambda(z)\notin\sigma(J)$, and let $J^{[m]}$ denote the Toeplitz-plus-finite-rank approximations of $J$ with principal resolvents $G^{[m]}$. Then $J^{[m]}\to J$ as $m\to\infty$, so by Theorem 5.14 there exists $M$ such that for all $m\geq M$, $\lambda(z)\notin\sigma(J^{[m]})$. For such $m$, both $G(\lambda(z))$ and $G^{[m]}(\lambda(z))$ are well defined, and using a well-known result on the difference of inverses (see for example, [2]), we have

[TABLE]

Theorem 5.13 shows that $\lim_{m\to\infty}c^{[m]}_{\mu}(z)/c^{[m]}(z)=c_{\mu}(z)/c(z)$. Therefore by Theorem 4.21 these limits are the same and we have equation (5.4). ∎

6 Computability aspects

In this section we discuss computability questions à la Ben-Artzi–Colbrook–Hansen–Nevanlinna–Seidel [4, 5, 28]. This involves an informal definition of the Solvability Complexity Index (SCI), a recent development that rigorously describes the extent to which various scientific computing problems can be solved. It is in contrast to classical computability theory à la Turing, in which problems are solvable exactly in finite time. In scientific computing we are often interested in problems which we can only approximate the solution in finite time, such that in an ideal situation this approximation can be made as accurate as desired. For example, the solution to a differential equation, the roots of a polynomial, or the spectrum of a linear operator.

Throughout this section we will consider only real number arithmetic, and the results do not necessarily apply to algorithms using floating point arithmetic.

The Solvability Complexity Index (SCI) has a rather lengthy definition, but in the end is quite intuitive [4].

Definition 6.1** (Computational problem).**

A computational problem is a 4-tuple, $\{\Xi,\Omega,\Lambda,\mathcal{M}\}$, where $\Omega$ is a set, called the input set, $\Lambda$ is a set of functions from $\Omega$ into the complex numbers, called the evaluation set, $\mathcal{M}$ is a metric space, and $\Xi:\Omega\to\mathcal{M}$ is the problem function.

Definition 6.2** (General Algorithm).**

Given a computational problem $\{\Xi,\Omega,\Lambda,\mathcal{M}\}$, a general algorithm is a function $\Gamma:\Omega\to\mathcal{M}$ such that for each $A\in\Omega$,

1. (i)

The action of $\Gamma$ on $A$ only depends on the set $\{f(A):f\in\Lambda_{\Gamma}(A)\}$ where $\Lambda_{\Gamma}(A)$ is a finite subset of $\Lambda$ 2. (ii)

For every $B\in\Omega$, $f(B)=f(A)$ for all $f\in\Lambda_{\Gamma}(A)$ implies $\Lambda_{\Gamma}(B)=\Lambda_{\Gamma}(A)$.

This definition of an algorithm is very general indeed. There are no assumptions on how $\Gamma$ computes its output; requirement $(i)$ ensures that it can only use a finite amount of information about its input, and requirement $(ii)$ ensures that the algorithm will only be affected by changes in the inputs which are actually measured. In short, $\Gamma$ depends only on, and is determined by, finitely many evaluable elements of each input.

Definition 6.3**.**

[Solvability Complexity Index] A computational problem function $\{\Xi,\Omega,\Lambda,\mathcal{M}\}$ has Solvability Complexity Index $k$ if $k$ is the smallest integer such that, for each $(n_{1},\ldots,n_{k})\in\mathbb{N}^{k}$ there exists a general algorithm $\Gamma_{n_{1},\ldots,n_{k}}:\Omega\to\mathcal{M}$, such that for all $A\in\Omega$,

[TABLE]

where the limit is taken in the metric on $\mathcal{M}$. In other words, the output of $\Gamma$ can be computed using a sequence of $k$ limits.

We require a metric space for the SCI.

Definition 6.4**.**

The Hausdorff metric for two compact subsets of the complex plane $A$ and $B$ is defined to be

[TABLE]

If a sequence of sets $A_{1},A_{2},A_{3},\ldots$ converges to $A$ in the Hausdorff metric, we write $A_{n}\xrightarrow{H}A$ as $n\to\infty$.

The computational problems considered for the remainder of this paper have $\Omega$ as a set of bounded self-adjoint operators on $\ell^{2}$, $\Lambda$ is the set of functions which simply return each individual element of the matrix representation of the operator, $\mathcal{M}$ is the set of subsets of $\mathbb{R}$ equipped with the Hausdorff metric, and $\Xi$ returns the spectrum of the operator.

Theorem 6.5** ([4]).**

The Solvability Complexity Index of the problem of computing the spectrum of a self-adjoint operator $A\in\mathcal{B}(\ell^{2})$ is equal to 2 with respect to the Hausdorff metric on $\mathbb{R}$. For compact operators and banded self-adjoint operators the SCI reduces to 1.

Theorem 6.5 implies that the SCI of computing the spectrum of bounded Jacobi operators in the Hausdorff metric is 1. In loose terms, the problem is solvable using only one limit of computable outputs. What more can we prove about the computability?

The results of Section 4 reduce the computation of the spectrum of a Toeplitz-plus-finite-rank Jacobi operator to finding the roots of a polynomial. From an uninformed position, one is lead to believe that polynomial rootfinding is a solved problem, with many standard approaches used every day. One common method is to use the QR algorithm to find the eigenvalues of the companion matrix for the polynomial. This can be done stably and efficiently in practice [3]. However, the QR algorithm is not necessarily convergent for non-normal matrices (companion matrices are normal if and only if they are unitary, which is exceptional). Fortunately, the SCI of polynomial rootfinding with respect to the Hausdorff metric in for subsets of $\mathbb{C}$ is 1, but if one requires the multiplicities of these roots then the SCI is not yet known [4].

A globally convergent polynomial rootfinding algorithm is given in [29]. For any degree $d$ polynomial the authors describe a procedure guaranteed to compute fewer than $1.11d(\log d)^{2}$ points in the complex plane, such that for each root of the polynomial, a Newton iteration starting from at least one of these points will converge to this root.

Let $\epsilon>0$. If a polynomial $p$ of degree $d$ has $r$ roots, how do we know when to stop so that we have $r$ points in the complex plane each within $\epsilon$ of a distinct root of $p$? This leads us to the concept of error control.

Definition 6.6**.**

[Error control] A function $\Gamma$ which takes inputs to elements in a metric space $\mathcal{M}$ is computable with error control if it has solvability complexity index 1, and for each $\epsilon$ we can compute $n$ to guarantee that

[TABLE]

In other words, the output of $\Gamma$ can be computed using a single limit, and an upper bound for the error committed by each $\Gamma_{n}$ is also computable.

Besides providing $\mathcal{O}(d(\log d)^{2})$ initial data for the Newton iteration (to find the complex roots of a degree $d$ polynomial), the authors of [29] discuss stopping criteria. In Section 9 of [29], it is noted therein that for Newton iterates $z_{1},z_{2},\ldots$, if $|z_{k}-z_{k-1}|<\epsilon/d$, then there exists a root $\xi$ of the polynomial in question such that $|z_{k}-\xi|<\epsilon$. It is then noted, however, that if there are multiple roots then it is in general impossible to compute their multiplicities with complete certainty. This is because the Newton iterates can pass arbitrarily close to a root to which this iterate does not, in the end, converge. Another consequence of this possibility is that roots could be missed out altogether because all of the iterates can be found to be close to a strict subset of the roots.

To salvage the situation, we give the following lemma, which adds some assumptions to the polynomial in question.

Lemma 6.7**.**

Let $p$ be a polynomial and $\Omega\subset\mathbb{C}$ an open set such that, a priori, the degree $d$ is known and it is known that there are $r$ distinct roots of $p$ in $\Omega$ and no roots on the boundary of $\Omega$. Then the roots of $p$ in $\Omega$ is computable with error control in the Hausdorff metric (see Definition 6.4 and Definition 6.6).

Proof.

Use Newton’s method with the $\mathcal{O}(d(\log d)^{2})$ complex initial data given in [29]. Using the stopping criteria in the discussion preceding this lemma, the algorithm at each iteration produces $\mathcal{O}(d(\log d)^{2})$ discs in the complex plane, within which all roots of $p$ must lie. To be clear, these discs have centres $z_{k}$ and radii $d\cdot|z_{k}-z_{k-1}|$. Let $R_{k}\subset\Omega$ denote the union of the discs which lie entirely inside $\Omega$ and have radius less than $\epsilon$ (the desired error). Note that this set may be empty if none of the discs are sufficiently small.

Because the Newton iterations are guaranteed to converge from these initial data, we must have eventually, for some sufficiently large $k$, that $R_{k}$ has $r$ connected components each with diameter less than $\epsilon$. Terminate when this verifiable condition has been fulfilled. ∎

Theorem 6.8**.**

Let $J=\Delta+F$ be a Toeplitz-plus-finite-rank Jacobi operator such that the rank of $F$ is known a priori. Then its point spectrum $\sigma_{p}(J)$ is computable with error control in the Hausdorff metric (see Definition 6.4 and Definition 6.6).

Remark 6.9*.*

Note that the full spectrum is simply $[-1,1]\cup\sigma_{p}(J)$.

Proof.

Suppose $F$ is zero outside the $n\times n$ principal submatrix. The value of $n$ can be computed given that we know the rank of $F$. Compute the principal $2n\times 2n$ submatrix of the connection coefficients matrix $C_{J\to\Delta}$ using formulae (3.3)–(3.7). The entries in the final column of this $2n\times 2n$ matrix give the coefficients of the Toeplitz symbol $c$, which is a degree $2n-1$ polynomial.

Decide if $\pm 1$ are roots by evaluating $p(\pm 1)$. Divide by the linear factors if necessary to obtain a polynomial $\tilde{p}$ such that $\tilde{p}(\pm 1)\neq 0$. Use Sturm’s Theorem to determine the number of roots of $\tilde{p}$ in $(-1,1)$, which we denote $r$ [39]. Since all roots in $\overline{\mathbb{D}}$ are real, there are $r$ roots of $\tilde{p}$ in the open unit disc $\mathbb{D}$ and none on the boundary.

By Lemma 6.7, the roots $z_{1},\ldots,z_{r}$ of this polynomial $c$ which lie in $(-1,1)$ can be computed with error control. By Theorem 4.22, for the point spectrum of $J$ we actually require $\lambda_{k}=\frac{1}{2}(z_{k}+z_{k}^{-1})$ to be computed with error control. Note that since $|\lambda_{k}|\leq\|J\|_{2}$ for each $k$, we have that $|z_{k}|\geq(1+2\|J\|_{2})^{-1}$. We should ensure that this holds for the computed roots $\hat{z}_{k}\in\mathbb{D}$ too. By the mean value theorem,

[TABLE]

Therefore it suffices to compute $\hat{z}_{k}$ such that $|z_{k}-\hat{z}_{k}|\leq\frac{\epsilon}{2}(1+\|F\|_{2})^{-1}(2+\|F\|_{2})^{-1}$, where $\epsilon$ is the desired error in the eigenvalues. ∎

The following Theorem shows that taking Toeplitz-plus-finite rank approximations of a Toeplitz-plus-compact Jacobi operator is sufficient for computing the spectrum with error control with respect to the Hausdorff metric.

Theorem 6.10**.**

Let $J=\Delta+K$ be a Toeplitz-plus-compact Jacobi operator. If for all $\epsilon>0$ an integer $m$ can be computed such that

[TABLE]

then the spectrum can be computed with error control in the Hausdorff metric.

Proof.

Let $\epsilon>0$. By the oracle assumed in the statement of the theorem, compute $m$ such that

[TABLE]

Now compute the point spectrum of the Toeplitz-plus-finite-rank approximation $J^{[m]}$ such that $d_{H}(\Sigma,\sigma(J^{[m]}))<\epsilon/2$, where $\Sigma$ denotes the computed set. Then, using Theorem 5.14, we have

[TABLE]

Here we used the fact that for a self-adjoint tridiagonal operator $A$,

[TABLE]

This completes the proof. ∎

An immediate question following Lemma 6.7 and Theorem 6.8 is why we have opted to use a Newton iteration in the complex plane instead of a purely real algorithm. We do this purely because Lemma 6.7 is an interesting point to make in and of itself with regards to the Solvability Complexity Index of polynomial rootfinding with error control. The key point is that while there exist algorithms to compute all of the roots of a polynomial (without multiplicity) in a single limit (i.e. with SCI equal to 1), one does not necessarily know when to stop the algorithm to achieve a desired error. Lemma 6.7 provides a basic condition on the polynomial to allow such control, which applies to this specific spectral problem.

7 Conclusions

In this paper we have proven new results about the relationship between the connection coefficients matrix between two different families of orthonormal polynomials, and the spectral theory of their associated Jacobi operators. We specialised the discussion to finite-rank perturbations of the free Jacobi operator and demonstrated explicit formulas for the principal resolvent and the spectral measure in terms of entries of the connection coefficients matrix. We showed that the results extend to trace class perturbations. Finally, we discussed computability aspects of the spectra of Toeplitz-plus-compact Jacobi operators. We showed that the spectrum of a Toeplitz-plus-compact Jacobi operator can be computed with error control, as long as the tail of the coefficients can be suitably estimated.

There are some immediate questions. Regarding regularity properties of the Radon-Nikodym derivative $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$ between the spectral measures $\nu$ and $\mu$ of Jacobi operators $D$ and $J$ respectively given in Propositions 3.6 and 3.9 and Corollary 3.10: can weaker regularity of $\frac{\mathrm{d}\nu}{\mathrm{d}\mu}$ be related to weak properties of $C=C_{J\to D}$? For example, the present authors conjecture that the Kullbeck–Leibler divergence,

[TABLE]

is finite if and only if the function of operators, $C^{T}C\log(C^{T}C)$ is well-defined as an operator mapping $\ell_{\mathcal{F}}\to\ell_{\mathcal{F}}^{\star}$. The reasoning comes from Lemma 3.8. Making such statements more precise for the case where $D=\Gamma$ or $D=\Delta$ (see equation (4.1)) could give greater insight into Sze̋go and quasi-Sze̋go asymptotics (respectively) for orthogonal polynomials [22, 14, 33].

Regarding computability: is there a theorem that covers the ground between Theorem 6.8 (for Toeplitz-plus-finite-rank Jacobi operators) and Theorem 6.10 (for Toeplitz-plus-compact Jacobi operators)? What can be said about the convergence of the continuous part of the spectral measure of a Toeplitz-plus-finite-rank truncations of a Toeplitz-plus-trace-class Jacob operator? Proposition 5.2 implies that this convergence is at least weak sense when tested against $f\in C_{b}(\mathbb{R})$.

The computability theorems in Section 6 all assume real arithmetic. What can be said about floating point arithmetic? Under what situations can the computation fail to give an unambiguously accurate solution? Answering this question is related to the mathematical problem of stability of the spectral measure under small perturbations of the Jacobi operator.

This paper also opens some broader avenues for future research. The connection coefficient matrix can be defined for any two Jacobi operators $J$ and $D$. It is natural to explore what structure $C_{J\to D}$ has when $D$ is a different reference operator to $\Delta$, and $J$ is a finite rank, trace class, or compact perturbation of $D$. For example, do perturbations of the Jacobi operator with periodic entries [12, 26] have structured connection coefficient matrices? Beyond periodic Jacobi operators, it would be interesting from the viewpoint of ergodic theory if we could facilitate the study and computation of almost-periodic Jacobi operators, such as the discrete almost-Mathieu operator [17]. Perturbations of the Jacobi operators for Laguerre polynomials and the Hermite polynomials could also be of interest, but challenges associated with the unboundedness of these operators could hamper progress [36]. Discrete Schrödinger operators with non-decaying potentials will also be of interest in this direction.

Spectra of banded self-adjoint operators may be accessible with these types of techniques too. Either using connection coefficient matrices between matrix orthogonal polynomials [13], or developing tridiagonalisation techniques are possible approaches, but the authors also consider this nothing more than conjecture at present. The multiplicity of the spectrum for operators with bandwidth greater than 1 appears to be a major challenge here. This becomes even more challenging for non-banded operators, such as Schrödinger operators on ${\mathbb{Z}}^{d}$ lattices.

Lower Hessenberg operators define polynomials orthogonal with respect to Sobolev inner products [24, pp. 40–43]. Therefore, we have two families of (Sobolev) orthogonal polynomials with which we may define connection coefficient matrices, as discussed in [27, p. 77]. Whether the connection coefficient matrices (which are still upper triangular) have structure which can be exploited for studying and computing the spectra of lower Hessenberg operators is yet to be studied.

Besides spectra of discrete operators defined on $\ell^{2}$, we conjecture that the results of this paper will also be applicable to continuous Schrödinger operators on $L^{2}(\mathbb{R})$, which are of the form $L_{V}[\phi](x)=-\phi^{\prime\prime}(x)+V(x)\phi(x)$ for a potential function $V:\mathbb{R}\to\mathbb{R}$. The reference operator is the negative Laplacian $L_{0}$ (which is the “free” Schrödinger operator). In this scenario, whereas the entries of a discrete connection coefficient matrix satisfy a discrete second order recurrence relation on $\mathbb{N}_{0}^{2}$ (see Lemma 3.2), the continuous analogue of the connection coefficient operator $C_{L_{V}\to L_{0}}$ is an integral operator whose (distributional) kernel satisfies a second order PDE on $\mathbb{R}^{2}$.

Appendix A Numerical results and the SpectralMeasures package

In this appendix we demonstrate some of the features of the SpectralMeasures.jl Julia package [53] that the authors have written to implement the ideas in the paper. This is part of the JuliaApproximation project, and builds on the package ApproxFun.jl [37]. ApproxFun is an extensive piece of software influenced by the Chebfun package [19] in Matlab, which can represent functions and operators [38, 37]. The code is subject to frequent changes and updates.

Given a Jacobi operator $J$ which is a finite-rank perturbation of the free Jacobi operator $\Delta$ with entries given by $\alpha_{k}=0$, $\beta_{k-1}=\frac{1}{2}$ for all $k\geq n$, SpectralMeasures.jl enables calculation of the following:

1. (i)

The connection coefficients matrix $C_{J\to\Delta}$: This is computed using the recurrences in equation (3.3)–(3.7). By Theorem 4.8, we only need to compute $n(n+1)$ entries of $C$ to have complete knowledge of all entries. In SpectralMeasures.jl, there is a type of operator called PertToeplitz, which allows such an operator to be stored and manipulated as if it were the full infinite-dimensional operator. 2. (ii)

The spectral measure $\mu(s)$: By Theorem 4.14, this measure has the form

[TABLE]

where $p_{C}$ is the polynomial given by the computable formula $p_{C}(s)=\sum_{k=0}^{2n-1}\langle C^{T}e_{k},C^{T}e_{0}\rangle U_{k}(s)$ and $r\leq n$. By Theorem 4.22, the numbers $\lambda_{k}$ are found by finding the distinct real roots $z_{k}$ of $c$ (the Toeplitz symbol of the Toeplitz part of $C$, which here is a polynomial of degree $2n-1$) in the interval $(-1,1)$. Also by Theorem 4.22, the weights $w_{k}$ can be computed using the formula

[TABLE] 3. (iii)

The principal resolvent $G(\lambda)$: For any $\lambda\in\mathbb{C}\setminus\sigma(J)$, by Theorem 4.12, this function can be defined by the formula

[TABLE]

where $p_{C}$ is as above and $p^{\mu}_{C}(\lambda)=\sum_{k=0}^{2n-1}\langle(C^{\mu})^{T}e_{k},C^{T}e_{0}\rangle U_{k}(\lambda)$. 4. (iv)

The mapped principal resolvent $G(\lambda(z))$: which is the principal resolvent mapped to $z$ in the unit disc by the Joukowski map $\lambda:z\to\frac{1}{2}(z+z^{-1})$. This is computed using the simple formula from Theorem 4.21,

[TABLE]

where $c$ and $c_{\mu}$ are the Toeplitz symbols of the Toeplitz parts of $C$ and $C^{\mu}$ respectively (these are polynomials of degree $2n-1$ and $2n-2$ respectively).

Consider a concrete example of a Toeplitz-plus-finite-rank Jacobi operator,

[TABLE]

The connection coefficients operator $C=C_{J\to\Delta}$ is

[TABLE]

As was noted above, we only need to explicitly compute $3\cdot 4=12$ entries and the rest are defined to be equal to the entry above and left one space. Because the perturbation is rational, the connection coefficients operator can be calculated exactly using the Rational type.

In Figure 2 we present plots of the spectral measure, the principal resolvent and the mapped principal resolvent. This format of figure is repeated for other Jacobi operators in the appendix.

The plot on left in Figure 2 is the spectral measure. There is a continuous part supported in the interval $[-1,1]$ and Dirac deltas are represented by vertical lines whose heights are precisely their weights. As noted in Theorem 6.8, it is possible to compute the eigenvalues of $J$ with guaranteed error control. Computations with guaranteed error control are made quite straightforward and flexible using the ValidatedNumerics.jl package [6], in which computations are conducted using interval arithmetic, and the desired solution is rigorously guaranteed to lie within the interval the algorithm gives the user [46]. Using this open source Julia package, we can compute the two eigenvalues for this operator to be $-1.1734766767874558$ and $1.5795946563898884$ with a guaranteed maximum error of $8.9\times 10^{-16}$. This can be replicated using the command validated_spectrum([.75;-.25;.5],[1;.75]) in SpectralMeasures.

The coloured plots in the middle and the right of Figure 2 are Wegert plots (sometimes called phase portraits) [54, 34]. For a function $f:\mathbb{C}\to\mathbb{C}$, a Wegert plot assigns a colour to every point $z\in\mathbb{C}$ by the argument of $f(z)$. Specifically, if $f(z)=re^{i\theta}$, then $\theta=0$ corresponds to the colour red, then cycles upwards through yellow, green, blue, purple as $\theta$ increases until at $\theta=2\pi$ it returns to red. This makes zeros and poles very easy to see, because around them the argument cycles through all the colours the same number of times as the degree of the root or pole. In these particular Wegert plots, we also plot lines of constant modulus as shadowed steps.

The middle plot in Figure 2 is the principal resolvent $G(\lambda)$, which always has a branch cut along the interval $[-1,1]$ and roots and poles along the real line. The poles correspond to Dirac delta measures in the spectral measure.

The third plot is the principal resolvent of $J$ mapped to the unit disc by the Joukowski map. Poles and roots of this resolvent in the unit disc correspond to those of the middle plot outside $[-1,1]$.

In Figure 3 we have plotted the spectral measure and principal resolvent of the Basic Perturbation 1 (see Examples 4.2, 4.17, 4.23) in which the top-left entry of the operator has been set to $\alpha/2$ for values $\alpha=0,0.15,0.35,0.5,0.75,1$. For the first four cases, the perturbation from the free Jacobi operator is small, and so the spectrum is purely continuous, which corresponds to no poles in the principal resolvent, and in the mapped resolvent there are only poles outside the unit disc. For the cases $\alpha=0.75$ and $1$, the Jacobi operator has a single isolated point of discrete spectrum. This is manifested as a Dirac delta in the spectral measure and a single pole in the principal resolvent.

In Figure 4 we have plotted the spectral measure and principal resolvent of Basic Perturbation 2 (see Examples 4.3, 4.18, 4.24) in which the $(0,1)$ and $(1,0)$ entries have been set to $\beta/2$ for values $\beta=0.5,0.707,0.85,1.0,1.2,1.5$. The effect is similar to that observed in Figure 3. For small perturbations the spectrum remains purely continuous, but for larger perturbations here two discrete eigenvalues emerge corresponding to Dirac deltas in the spectral measure and poles in the principal resolvent.

In Figure 5 we have plotted a sequence of approximations to the Jacobi operator for the Legendre polynomials, which has entries $\alpha_{k}=0$ for $k=0,1,2,\ldots$ and

[TABLE]

This is a Toeplitz-plus-trace-class Jacobi operator because $\beta_{k}=\frac{1}{2}+\mathcal{O}(k^{-2})$, and by taking Toeplitz-plus-finite-rank approximations $J^{[n]}$ as in equation (5.1), we can compute approximations to the spectral measure and principal resolvent. Figure 5 depicts the spectral measure and the principal resolvent for the Toeplitz-plus-finite-rank Jacobi operators $J^{[n]}$ for the values $n=1,2,3,10,30,100$. For the spectral measures, we see that there is no discrete part for any $n$, and as $n$ increases, the spectral measure converges to the scaled Lebesgue measure $\frac{1}{2}\mathrm{d}s$ restricted to $[-1,1]$. The convergence is at least weak by Proposition 5.2, but it would be interesting (as mentioned in the conclusions) to determine if there is a stronger form of convergence at play due to the perturbation of $\Delta$ lying in the space of trace class operators. There is a Gibbs-like effect occurring at the boundaries, which suggests that this convergence occur pointwise everywhere up to the boundary of $[-1,1]$. For the principal resolvents, the middle plots do not vary greatly, as the difference between the functions in the complex plane is not major. However, in right plots, there are hidden pole-root pairs in the resolvent lying outside the unit disc which coalesce around the unit disc and form a barrier. The meaning of this barrier is unknown to the authors.

Figures 6 and 7 demonstrate similar features to Figure 5, except that the polynomials sequences they correspond to are the ultraspherical polynomials with parameter $\gamma=0.6$ (so that the spectral measure is proportional to $(1-s^{2})^{1.1}$) and the Jacobi polynomials with parameters $(\alpha,\beta)=(0.4,1.9)$ (so that the spectral measure is proportional to $(1-s)^{0.4}(1+s)^{1.9}$). Similar barriers of pole-root pairs outside the unit disc occur for these examples as well.

Figure 8 presents a Toeplitz-plus-trace-class Jacobi operator with pseudo-randomly generated entries. With a random vector $\mathbf{r}$ containing entries uniformly distributed in the interval $[0,1)$, the following entries were used

[TABLE]

Then the Toeplitz-plus-finite-rank approximations $J^{[n]}$ (see equation (5.1)) of this operator were taken for values $n=1,2,3,10,50,100$. Since the off-diagonal elements are constant, this is a scaled and shifted version of a discrete Schrödinger operator with a random, decaying potential.