# Fast non-Hermitian Toeplitz eigenvalue computations, joining matrix-less
algorithms and FDE approximation matrices

M. Bogoya<EMAIL_ADDRESS>Universidad del Valle, Cali,
Colombia.
Departamento de Matemáticas. S.M. Grudsky<EMAIL_ADDRESS>CINVESTAV,
Mexico, Mexico D.F.
Departamento de Matemáticas.  S. Serra–Capizzano
<EMAIL_ADDRESS>University of Insubria, Como, Italy.
Dipartimento di Scienza e Alta Tecnologia.

###### Abstract

The present work is devoted to the eigenvalue asymptotic expansion of the
Toeplitz matrix $T_{n}(a)$ whose generating function $a$ is complex valued and
has a power singularity at one point. As a consequence, $T_{n}(a)$ is non-
Hermitian and we know that the eigenvalue computation is a non-trivial task in
the non-Hermitian setting for large sizes. We follow the work of Bogoya,
Böttcher, Grudsky, and Maximenko and deduce a complete asymptotic expansion
for the eigenvalues. After that, we apply matrix-less algorithms, in the
spirit of the work by Ekström, Furci, Garoni, Serra-Capizzano et al, for
computing those eigenvalues. Since the inner and extreme eigenvalues have
different asymptotic behaviors, we worked on them independently, and combined
the results to produce a high precision global numerical and matrix-less
algorithm.

The numerical results are very precise and the computational cost of the
proposed algorithms is independent of the size of the considered matrices for
each eigenvalue, which implies a linear cost when all the spectrum is
computed. From the viewpoint of real world applications, we emphasize that the
matrix class under consideration includes the matrices stemming from the
numerical approximation of fractional diffusion equations. In the final
conclusion section a concise discussion on the matter and few open problems
are presented.
MSC Classes: Primary 15B05, 65F15, 47B35. Secondary 15A18, 47A38.

##  1 Introduction

Let $\mathbb{T}$ be the complex unit circle. For a complex-valued function
$a\in L^{1}(\mathbb{T})$ let $a_{k}$ be its $k$th Fourier coefficient, that is

$a_{k}\equiv\frac{1}{2\pi}\int_{0}^{2\pi}a(\hskip 0.6pt\mathrm{e}^{\hskip
0.2pt\mathrm{i}\theta})\hskip 0.6pt\mathrm{e}^{-\hskip
0.2pt\mathrm{i}k\theta}\operatorname{d\\!}\theta,\quad\ k\in\mathbb{Z},\
\hskip 0.2pt\mathrm{i}^{2}=-1.$

The $n\times n$ Toeplitz matrix with generating function $a$ is given by
$T_{n}(a)\equiv(a_{j-k})_{j,k=0}^{n-1}$, and is characterized by having the
fixed constant $a_{k}$ along the $k$th diagonal, parallel to the main one. In
formulae, we have

$T_{n}(a)=\begin{bmatrix}a_{0}&a_{-1}&a_{-2}&\cdots&\cdots&a_{-(n-1)}\\\
a_{1}&\ddots&\ddots&\ddots&&\vdots\\\
a_{2}&\ddots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&a_{-2}\\\
\vdots&&\ddots&\ddots&\ddots&a_{-1}\\\
a_{n-1}&\cdots&\cdots&a_{2}&a_{1}&a_{0}\end{bmatrix}.$

We will follow the notation of [1]. Let $H^{\infty}$ be the Hardy space of
(boundary values of) bounded analytic functions in the complex unit disk
$\mathbb{D}$, $\mathbb{T}=\partial\mathbb{D}$. For a continuous generating
function $a$, we denote by

$\operatorname{wind}_{\lambda}(a)\equiv\frac{1}{2\pi\hskip
0.2pt\mathrm{i}}\oint_{\mathcal{R}(a)}\frac{\operatorname{d\\!}\zeta}{\zeta-\lambda},$

the winding number of $a$ in connection with a point
$\lambda\in\mathbb{C}\setminus\mathcal{R}(a)$ where $\mathcal{R}(a)$ is the
range of $a$. Intuitively, $\operatorname{wind}_{\lambda}(a)$ is the number of
times that the range of $a$ travels around the point $\lambda$ in
counterclockwise sense. We also denote by $\mathcal{D}(a)$ the set of all
$\lambda\in\mathbb{C}$ for which $\operatorname{wind}_{\lambda}(a)\neq 0$. In
plane words $\mathcal{D}(a)$ is the collection of all the points trapped by
$\mathcal{R}(a)$. In the first part of the current work, we will study the
eigenvalues of $T_{n}(a)$ for functions $a(t)=\frac{1}{t}(1-t)^{\alpha}f(t)$
satisfying the following properties:

  1. (i)

The constant $\alpha$ belongs to the interval $(0,1)$ and $f\in
C^{\infty}(\mathbb{T})\cap H^{\infty}$.

  2. (ii)

$f(z)\neq 0$ for $z\in\overline{\mathbb{D}}$, and $f$ has an analytic
extension to an open neighborhood $W$ of $\mathbb{T}\setminus\\{1\\}$ not
containing the point $1$.

  3. (iii)

The range of $a$ is a Jordan curve in $\mathbb{C}$ which is the boundary of a
connected region in the complex plane.

Strictly speaking, the complex-valued function $(1-z)^{\alpha}$ is multi-
valued, hence we take the argument specified by
$\arg(1-z)^{\alpha}\in(-\alpha\pi,\alpha\pi]$ when $\arg(1-z)\in(-\pi,\pi]$.
Property (ii) is necessary because, in that case, $b$ belongs to the Hardy
space $H^{\infty}$ and the nice relation
$T_{n}^{-1}(b)=T_{n}\big{(}\frac{1}{b}\big{)}$ holds, which is used for
expressing the determinant $D_{n}(a)$ of $T_{n}(a)$ as a Fourier integral.
Additionally, the condition $\alpha\in(0,1)$ implies that $a$ has a power
singularity at $t=1$, i.e. no continuous derivatives there.

At this point it is worth noticing that the considered set of functions
includes the generating functions of a class of Hessenberg Toeplitz matrices
stemming from the numerical approximation of certain fractional diffusion
equations (FDEs); see [5, 6, 13]. In the latter context, the parameter
$\alpha$ belongs to the open interval $(0,2)$. We remind that the quoted FDEs
with noninteger parameter $\alpha\in(0,2)$ have gained a tremendous attention
in real world applications, because they model anomalous diffusion processes
arising in Biology, Physics, etc. (see again [13] and references therein).

Property (ii) also guaranties that the only point of $\mathbb{T}$ where $a$ is
not analytic is $t=1$, and finally, property (iii) amounts in avoiding that
$\mathcal{R}(a)$ shows the presence of loops. The resulting Toeplitz matrices
$T_{n}(a)$ are lower Hessenberg, meaning that they can be obtained by adding
one upper diagonal to a lower triangular matrix.

We emphasize that this knowledge is important when analyzing fast methods for
the solution of large linear systems arising from the numerical approximation
of the considered FDEs, via preconditioned Krylov methods, multigrid
algorithms, and multi-iterative solvers (see e.g. [12] and references
therein).

In [1] an asymptotic expansion was deduced for the eigenvalues bounded away
from zero, which the authors named inner, while in [4] the extreme (closest to
zero) eigenvalues were studied. With a different technique, the authors were
able to obtain an asymptotic expression for the eigenvalues located near to
the singularity point. The result is reminiscent of the works by S. Parter
[18, 19] and H. Widom [24], for the real-valued case; see also [9, 20, 21, 23]
for a more complete picture on the subject, including newer matrix theoretic
tools.

The case $\alpha\in(1,2)$ or more generally,
$\alpha\in\mathbb{R}_{+}\setminus\mathbb{Z}$, can be managed analogously for
the inner eigenvalues, while the extreme eigenvalues remain a challenge that
will be the subject of a future investigation. For this reason we restrict
ourselves to the case $\alpha\in(0,1)$.

Regarding computational costs, the best option relies in employing
extrapolation techniques based on the asymptotic expansion, because the latter
leads to a direct eigenvalue approximation which is independent of the matrix
entries or its size. Indeed the computation can be done in two different ways:
the first one is to use the analytic expression of the coefficient functions,
if available; the second one is to employ a numerical and matrix-less
algorithm to approximate them. After the respective precomputing phases, in
both cases, the calculation will need only a fraction of a second in a
standard computer. The standard option is to use any commercial eigensolver
which are limited by the matrix size and its conditioning. To put it in
numbers, the eigenvalue calculation of a $4000$ size moderately conditioned
matrix with 50 precision digits, can be done in a week by a standard computer.
If the matrix size is $10000$, only a supercomputer can handle the same task,
but if the matrix size is $10^{10}$ (a size appearing in several applications)
no computer can do it.

Taking inspiration from the literature, in the present work we split our study
into the mentioned cases. The article is organized as follows: in Section 2 we
deduce a complete asymptotic expansion for the inner eigenvalues of $T_{n}(a)$
which served as the theoretical ground for the subsequent matrix-less
numerical algorithm. In Section 3 we use a previous theoretical result to
conjecture an expansion which was the base for the respective matrix-less
numerical algorithm. In Section 4 we join the two processes in order to design
a global algorithm approximating all the eigenvalues with high accuracy: in
this section we also report and discuss our numerical experiments. Finally, in
Section 5 we draw some future works.

##  2 Inner eigenvalues

In the paper [1], the authors employed the Baxter–Schmidt formula (see [10,
§2.3]) for Toeplitz determinants to express $D_{n}(a-\lambda)$ as a Fourier
integral. In particular they also showed that if $W_{0}$ is a small open
neighborhood of zero in $\mathbb{C}$, for every point
$\lambda\in\mathcal{D}(a)\cap(a(W)\setminus W_{0})$, there exists a unique
point $t_{\lambda}\in\mathbb{C}$, with $|t_{\lambda}|>1$ such that
$a(t_{\lambda})=\lambda$. See Figure 1.

In order to attain an asymptotic expansion for the eigenvalues of $T_{n}(a)$,
we are based on [1] where the authors obtained an asymptotic expression
involving only two terms, our aim is to obtain an expression with any number
of terms. The following is our main result.

Figure 1: The range and neighborhoods of the generating function
$a(t)=\frac{1}{t}(1-t)^{\alpha}$ with $\alpha=\frac{3}{4}$. Left: For $n=12$
the red dots are the roots of unity $\omega_{n}^{j}$, the cyan dots are the
sequence $(t_{j,n})$ described in Theorem 2.1, that is
$\lambda_{j}(T_{n}(a))=a(t_{j,n})$, $j=0,\ldots,11$. The gray shaded region is
the neighborhood $W\setminus\overline{\mathbb{D}}$ where $a$ has an analytic
extension and where the sequence $t_{j,n}$ is located. Right: For $n=12$ the
cyan dots are the eigenvalues $\lambda_{j}(T_{12}(a))=\lambda_{j,12}$,
$j=1,\ldots,12$, the gray shaded region is the neighborhood
$a(W\setminus\overline{\mathbb{D}})$ where the inner eigenvalues are located.
The light red shaded region is the neighborhood $W_{0}=B(0,\varepsilon)$,
where the extreme eigenvalues are located.

######  Theorem 2.1.

Let $a(t)=\frac{1}{t}b(t)$ be a function satisfying the properties (i)–(iii).
Then for every small open neighborhood $W_{0}$ of zero in $\mathbb{C}$ and
every index $j$ with $a(\omega_{n}^{j})\notin W_{0}$, $t_{j,n}\equiv
t_{\lambda_{j,n}}$ admits the following asymptotic expansion

$\displaystyle t_{j,n}$ $\displaystyle\sim$
$\displaystyle\omega_{n}^{j}n^{(\alpha+1)\frac{1}{n}}\bigg{\\{}1+\sum_{r=0}^{\infty}\sum_{m=1}^{\infty}\sum_{\ell=0}^{m-1}p_{r,m,\ell}(\omega_{n}^{j})\frac{\log^{\ell}(n)}{n^{\alpha
r+m}}\bigg{\\}}$ $\displaystyle=$
$\displaystyle\omega_{n}^{j}\bigg{\\{}1+\sum_{r=0}^{\infty}\sum_{m=1}^{\infty}\sum_{\ell=0}^{m-1}q_{r,m,\ell}(\omega_{n}^{j})\frac{\log^{\ell}(n)}{n^{\alpha
r+m}}\bigg{\\}},$

where

  * •

$\lambda_{j}(T_{n}(a))=a(t_{j,n})$ and $\omega_{n}\equiv\exp(-\frac{2\pi\hskip
0.2pt\mathrm{i}}{n})$;

  * •

the coefficients $p_{r,m,\ell}$ and $q_{r,m,\ell}$ are functions in
$C(\mathbb{T})$ depending only on the generating function $a$.

To help the visualization of the expansion in the previous theorem, take
$h\equiv\frac{1}{n}$ and arrange the involved powers of $h$ with the following
array:

$h$ |  |  | $h^{\alpha+1}$ |  |  | $\cdots$
---|---|---|---|---|---|---
$h^{2}\log(n)$ | $h^{2}$ |  | $h^{\alpha+2}\log(n)$ | $h^{\alpha+2}$ |  | $\cdots$
$h^{3}\log^{2}(n)$ | $h^{3}\log(n)$ | $h^{3}$ | $h^{\alpha+3}\log^{2}(n)$ | $h^{\alpha+3}\log(n)$ | $h^{\alpha+3}$ | $\cdots$
$\vdots$ |  |  | $\vdots$ |  |  | $\ddots$

Hence the terms can be ordered, but the ordering depends on $\alpha$. In [1,
Th. 1.2] the coefficient $p_{0,1,0}$ was exactly calculated as

$p_{0,1,0}(z)=\,\log\Big{(}\frac{a^{2}(z)}{c\,a^{\prime}(z)z^{2}}\Big{)}\quad\mbox{with}\quad
c=\frac{1}{\pi}f(1)\Gamma(\alpha+1)\sin(\alpha\pi),$

where $\Gamma$ is the well-known Gamma function. In this section, instead of
the exact calculation of the remaining coefficients $p_{r,m,\ell}$ or
$q_{r,m,\ell}$, we are interested in justifying the whole expansion and then,
use it for proposing a matrix-less numerical algorithm, in the same fashion as
in [3, 8, 7, 14, 15, 16]. Indeed, we remind that the quoted procedures are all
essentially of extrapolation type and hence the formal (or numerical)
expression of the coefficients is not required, likewise the extrapolation
techniques for the computation of definite integrals in the sense of Romberg
[22, §3.4].

In asymptotic analysis we say that a function $f$ admits the asymptotic
expansion $\sum_{k=0}^{\infty}a_{k}(x-x_{0})^{k}$ at $x_{0}$, denoted as

$f(x)\sim\sum_{k=0}^{\infty}a_{k}(x-x_{0})^{k},$

if we have $f(x)=\sum_{k=0}^{m-1}a_{k}(x-x_{0})^{k}+O(|x-x_{0}|^{m})$ for
every $m\in\mathbb{Z}_{+}$. For proving Theorem 2.1 we will need some
preliminary findings. We start with the following well-known result, which is,
for example in [17, p. 97]. We present it here for reader convenience and
without a proof.

######  Theorem 2.2.

Let $\gamma>0$, $\delta>0$, $u\colon[0,\delta]\to\mathbb{C}$ be a $C^{\infty}$
function with $u^{(m)}(\delta)=0$ for all $m\geqslant 0$. Then, as
$n\to\infty$

$\int_{0}^{\delta}\theta^{\gamma-1}u(\theta)\hskip 0.6pt\mathrm{e}^{\hskip
0.2pt\mathrm{i}n\theta}\operatorname{d\\!}\theta\sim\sum_{m=0}^{\infty}\frac{u^{(m)}(0)\Gamma(m+\gamma)\hskip
0.6pt\mathrm{e}^{\hskip
0.2pt\mathrm{i}\frac{\pi}{2}(m+\gamma)}}{m!\,n^{m+\gamma}}.$

For a generating function $a(t)=\frac{1}{t}b(t)$ satisfying the properties
(i)–(iii), we can write

$b(\hskip 0.6pt\mathrm{e}^{\hskip 0.2pt\mathrm{i}\theta})=(1-\hskip
0.6pt\mathrm{e}^{\hskip 0.2pt\mathrm{i}\theta})^{\alpha}f(\hskip
0.6pt\mathrm{e}^{\hskip 0.2pt\mathrm{i}\theta})=\hskip
0.6pt\mathrm{e}^{-\frac{\hskip
0.2pt\mathrm{i}}{2}\alpha\pi}\theta^{\alpha}f_{0}(\theta)f(\hskip
0.6pt\mathrm{e}^{\hskip 0.2pt\mathrm{i}\theta}),$

where $f_{0}$ belongs to the class $C^{\infty}[-\delta,\delta]$. Hence $b$ has
a zero of order $\alpha$ at $\theta=0$. Let $\delta$ be a small positive
constant, $\phi\colon\mathbb{R}\to[0,1]$ a $C^{\infty}$ function with support
in $[-\delta,\delta]$, and $\phi|_{[-\varepsilon,\varepsilon]}\equiv 1$ for
some $0<\varepsilon<\delta$. For each $r\in\mathbb{Z}_{+}$ consider the
auxiliary function

$u_{r}(\theta)\equiv f_{0}^{r}(-\theta)f^{r}(\hskip 0.6pt\mathrm{e}^{-\hskip
0.2pt\mathrm{i}\theta})\phi(\theta)\hskip 0.6pt\mathrm{e}^{\hskip
0.2pt\mathrm{i}\theta(r+1)},$ (2.1)

and let $D_{n}(a)$ be the determinant of the Toeplitz matrix $T_{n}(a)$. The
following result gives us an expansion for the determinant of
$T_{n}(a-\lambda)$, which will be the starting point for obtaining the
expansion of the eigenvalues.

######  Theorem 2.3.

Let $a(t)=\frac{1}{t}b(t)$ be a function satisfying the properties (i)–(iii).
Then the determinant of $T_{n}(a-\lambda)$ is given by

$D_{n}(a-\lambda)=(-1)^{n}b_{0}^{n+1}\beta_{n}(\lambda),$

where $\beta_{n}$ admits the asymptotic expansion

$\beta_{n}(\lambda)\sim\frac{-1}{t_{\lambda}^{n+2}a^{\prime}(t_{\lambda})}-\frac{1}{\pi}\sum_{r,m=1}^{\infty}\frac{\Gamma(\alpha
r+m)\mathfrak{\,Re}\\{u_{r}^{(m-1)}(0)\hskip 0.6pt\mathrm{e}^{\frac{\hskip
0.2pt\mathrm{i}}{2}\pi(2\alpha r+m)}\\}}{(m-1)!a^{r+1}(t_{\lambda})n^{\alpha
r+m}},$

$b_{0}\neq 0$ is the zeroth Fourier coefficient of $b$, and $u_{r}$ is given
by (2.1).

###### Proof.

Recall that $\mathcal{D}(a)$ is the set of all complex points $z$ with
$\operatorname{wind}_{\lambda}(z)\neq 0$. From [1, Lm. 4.1 and Lm. 4.2], for
any $\lambda\in\mathcal{D}(a)$, we reach
$D_{n}(a-\lambda)=(-1)^{n}b_{0}^{n+1}\beta_{n}(\lambda)$, where $b_{0}\neq 0$
is the zeroth Fourier coefficient of $b$ and $\beta_{n}$ is a linear
combination of Fourier-type integrals, admitting the representation

$\beta_{n}(\lambda)=\frac{-1}{t_{\lambda}^{n+2}a^{\prime}(t_{\lambda})}-\frac{1}{2\pi}\sum_{r=1}^{\infty}\frac{I(r,n)}{a^{r+1}(t_{\lambda})}+O(n^{-\infty}),$
(2.2)

the order term is uniform in $\lambda$ and

$I(r,n)\equiv\int_{-\delta}^{\delta}b^{r}(\hskip 0.6pt\mathrm{e}^{\hskip
0.2pt\mathrm{i}\theta})\phi(\theta)\hskip 0.6pt\mathrm{e}^{-\hskip
0.2pt\mathrm{i}\theta(r+1)}\hskip 0.6pt\mathrm{e}^{-\hskip
0.2pt\mathrm{i}n\theta}\operatorname{d\\!}\theta.$

We will expand the integral $I$ with the aid of Theorem 2.2. Note that

$\displaystyle I(r,n)$ $\displaystyle=$ $\displaystyle\hskip
0.6pt\mathrm{e}^{-\frac{\hskip 0.2pt\mathrm{i}}{2}\alpha
r\pi}\int_{-\delta}^{\delta}\theta^{\alpha r}u_{r}(-\theta)\hskip
0.6pt\mathrm{e}^{-\hskip 0.2pt\mathrm{i}n\theta}\operatorname{d\\!}\theta$
$\displaystyle=$ $\displaystyle\hskip 0.6pt\mathrm{e}^{-\frac{\hskip
0.2pt\mathrm{i}}{2}\alpha
r\pi}\bigg{\\{}\overline{\int_{0}^{\delta}\theta^{\alpha
r}\overline{u_{r}(-\theta)}\hskip 0.6pt\mathrm{e}^{\hskip
0.2pt\mathrm{i}n\theta}\operatorname{d\\!}\theta}+\hskip
0.6pt\mathrm{e}^{\hskip 0.2pt\mathrm{i}\alpha\pi
r}\int_{0}^{\delta}\theta^{\alpha r}u_{r}(\theta)\hskip
0.6pt\mathrm{e}^{\hskip
0.2pt\mathrm{i}n\theta}\operatorname{d\\!}\theta\bigg{\\}},$

where $u_{r}$ is given by (2.1). Using the relation
$\overline{u_{r}(\theta)}=u_{r}(-\theta)$ and applying Theorem 2.2 twice, we
obtain

$\displaystyle I(r,n)$ $\displaystyle\sim$ $\displaystyle\hskip
0.6pt\mathrm{e}^{-\frac{\hskip 0.2pt\mathrm{i}}{2}\alpha
r\pi}\bigg{\\{}\overline{\sum_{m=1}^{\infty}\frac{\Gamma(\alpha
r+m)}{(m-1)!}u_{r}^{(m-1)}(0)\hskip 0.6pt\mathrm{e}^{\frac{\hskip
0.2pt\mathrm{i}}{2}\pi(\alpha r+m)}n^{-\alpha r-m}}$ $\displaystyle+\hskip
0.6pt\mathrm{e}^{\hskip 0.2pt\mathrm{i}\alpha\pi
r}\sum_{m=1}^{\infty}\frac{\Gamma(\alpha r+m)}{(m-1)!}u_{r}^{(m-1)}(0)\hskip
0.6pt\mathrm{e}^{\frac{\hskip 0.2pt\mathrm{i}}{2}\pi(\alpha r+m)}n^{-\alpha
r-m}\bigg{\\}}$ $\displaystyle=$ $\displaystyle
2\sum_{m=1}^{\infty}\frac{\Gamma(\alpha
r+m)}{(m-1)!}\mathfrak{\,Re}\\{u_{r}^{(m-1)}(0)\hskip
0.6pt\mathrm{e}^{\frac{\hskip 0.2pt\mathrm{i}}{2}\pi(2\alpha
r+m)}\\}n^{-\alpha r-m}.$

Finally, replacing in (2.2), we attain

$\beta_{n}(\lambda)\sim\frac{-1}{t_{\lambda}^{n+2}a^{\prime}(t_{\lambda})}-\frac{1}{\pi}\sum_{r,m=1}^{\infty}\frac{\Gamma(\alpha
r+m)}{(m-1)!a^{r+1}(t_{\lambda})}\mathfrak{\,Re}\\{u_{r}^{(m-1)}(0)\hskip
0.6pt\mathrm{e}^{\frac{\hskip 0.2pt\mathrm{i}}{2}\pi(2\alpha
r+m)}\\}n^{-\alpha r-m},$

which proves the assertion of the theorem. ∎

The next step consists in solving the relation $D_{n}(a-\lambda)=0$ for
$t_{\lambda}$, which, with the aid of Theorem 2.3, leads to

$t_{\lambda}\sim\bigg{\\{}\frac{t_{\lambda}^{2}a^{\prime}(t_{\lambda})}{a^{2}(t_{\lambda})}\sum_{r,m=0}^{\infty}\frac{\varrho(r,m)}{a^{r}(t_{\lambda})n^{\alpha(r+1)+(m+1)}}\bigg{\\}}^{-\frac{1}{n}},$

where

$\varrho(r,m)\equiv-\frac{\Gamma(\alpha(r+1)+(m+1))}{\pi
m!}\mathfrak{\,Re}\\{u_{r}^{(m)}(0)\hskip 0.6pt\mathrm{e}^{\frac{\hskip
0.2pt\mathrm{i}}{2}\pi(2\alpha(r+1)+(m+1))}\\}.$

In order to expand the above expression and determine the powers of $n$
involved, we note that
$\varrho(0,0)=\frac{1}{\pi}\Gamma(\alpha+1)\sin(\alpha\pi)f(1)\neq 0$ and
consequently we write

$\displaystyle t_{\lambda}$ $\displaystyle\sim$
$\displaystyle\omega_{n}^{j}\bigg{\\{}\frac{\varrho(0,0)t_{\lambda}^{2}a^{\prime}(t_{\lambda})}{a^{2}(t_{\lambda})n^{\alpha+1}}\bigg{\\}}^{-\frac{1}{n}}\bigg{\\{}1+\frac{1}{\varrho(0,0)}\sum_{\begin{subarray}{c}r,m=0\\\
(r,m)\neq(0,0)\end{subarray}}^{\infty}\frac{\varrho(r,m)}{a^{r}(t_{\lambda})n^{\alpha
r+m}}\bigg{\\}}^{-\frac{1}{n}}$ (2.3) $\displaystyle=$
$\displaystyle\omega_{n}^{j}n^{(\alpha+1)\frac{1}{n}}\bigg{\\{}\frac{a^{2}(t_{\lambda})}{\varrho(0,0)t_{\lambda}^{2}a^{\prime}(t_{\lambda})}\bigg{\\}}^{\frac{1}{n}}\bigg{\\{}1+\frac{1}{\varrho(0,0)}\sum_{\begin{subarray}{c}r,m=0\\\
(r,m)\neq(0,0)\end{subarray}}^{\infty}\frac{\varrho(r,m)}{a^{r}(t_{\lambda})n^{\alpha
r+m}}\bigg{\\}}^{-\frac{1}{n}},$

where $j=0,\ldots,n-1$ and $\omega_{n}=\exp(-\frac{2\pi\hskip
0.2pt\mathrm{i}}{n})$. Now we want to expand the previous expression until the
outer exponents $\pm\frac{1}{n}$ disappear. The main aim here is to determine
the powers of $n$ involved instead of the value of the respective exponents.
One possible strategy for this aim is as follows. Take $h=\frac{1}{n}$. The
expansion of the first term involves $h$ with all natural exponents, i.e.
$h,h^{2},h^{3},h^{4},\ldots$, while the expansion of the second term involves
all the exponents in the array

$M\equiv\begin{bmatrix}&h^{\alpha}&h^{2\alpha}&h^{3\alpha}&\cdots\\\
h&h^{\alpha+1}&h^{2\alpha+1}&h^{3\alpha+1}&\cdots\\\
h^{2}&h^{\alpha+2}&h^{2\alpha+2}&h^{3\alpha+2}&\cdots\\\
h^{3}&h^{\alpha+3}&h^{2\alpha+3}&h^{3\alpha+3}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{bmatrix},$

combined as $1+hM+h(h+1)M^{2}+h(h+1)(h+2)M^{3}+\ldots$. Note that the entire
powers of the array $M$ will involve powers of $h$ following the next pattern

$M^{2}=\begin{bmatrix}&&h^{2\alpha}&h^{3\alpha}&\cdots\\\
&h^{\alpha+1}&h^{2\alpha+1}&h^{3\alpha+1}&\cdots\\\
h^{2}&h^{\alpha+2}&h^{2\alpha+2}&h^{3\alpha+2}&\cdots\\\
h^{3}&h^{\alpha+3}&h^{2\alpha+3}&h^{3\alpha+3}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{bmatrix},\
M^{3}=\begin{bmatrix}&&&h^{3\alpha}&\cdots\\\
&&h^{2\alpha+1}&h^{3\alpha+1}&\cdots\\\
&h^{\alpha+2}&h^{2\alpha+2}&h^{3\alpha+2}&\cdots\\\
h^{3}&h^{\alpha+3}&h^{2\alpha+3}&h^{3\alpha+3}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{bmatrix}.$

Hence, after expanding, multiplying, and ordering, we seize that the second
factor in (2.3) involves the constant $1$ plus terms involving $h$ with all
the powers in the set $\\{\alpha r+m\colon r,m\in\mathbb{Z},\,r\geqslant
0,\,m\geqslant 1\\}$. Finally, combining the expansions of the two factors of
(2.3) we obtain

$t_{\lambda}\sim\omega_{n}^{j}n^{(\alpha+1)\frac{1}{n}}\bigg{\\{}1+\sum_{\begin{subarray}{c}r\geqslant
0\\\ m\geqslant 1\end{subarray}}\frac{\xi_{r,m}(t_{\lambda})}{n^{\alpha
r+m}}\bigg{\\}},$ (2.4)

where $\xi_{r,m}$ are functions in
$C^{\infty}(W\setminus\overline{\mathbb{D}})$, that can be determined
explicitly. The previous relation is an implicit equation for $t_{\lambda}$.
In [1, Th. 5.1] it was proved that there exists a collection of mutually
disjoint neighborhoods of each $\omega_{n}^{j}$, $j=0,\ldots,n-1$, containing
exactly one solution of (2.4). We are now ready to prove our main result.

###### Proof of Theorem 2.1.

Let $t_{\lambda}$ be a solution of (2.4). Since $a(t_{\lambda})=\lambda$ must
be an eigenvalue of $T_{n}(a)$, we conclude that there are exactly $n$
solutions and we denote them as $t_{j,n}\equiv t_{\lambda_{j,n}}$,
$j=1,\ldots,n$. Now we iterate over (2.4). After expanding the term
$n^{(\alpha+1)\frac{1}{n}}$ and multiplying, we reach

$t_{\lambda}\sim\omega_{n}^{j}\bigg{\\{}1+\sum_{\ell=1}^{\infty}\zeta^{(1)}_{l}(t_{\lambda})\frac{\log^{\ell}(n)}{n^{\ell}}+\sum_{r=0}^{\infty}\sum_{m=1}^{\infty}\sum_{\ell=0}^{m-1}\zeta^{(2)}_{r,m,l}(t_{\lambda})\frac{\log^{\ell}(n)}{n^{\alpha
r+m}}\bigg{\\}},$ (2.5)

where the functions $\zeta_{l}^{(1)}$ and $\zeta_{r,m.l}^{(2)}$ belong to
$C^{\infty}(W\setminus\overline{\mathbb{D}})$ and can be determined
explicitly. Since the expanded terms in (2.5) have order $\frac{\log(n)}{n}$,
our first iteration is
$t_{\lambda}^{(1)}=\omega_{n}^{j}+O(\frac{\log(n)}{n})$. The second iteration
is

$\displaystyle t_{\lambda}^{(2)}$ $\displaystyle\sim$
$\displaystyle\omega_{n}^{j}n^{(\alpha+1)\frac{1}{n}}\bigg{\\{}1+\sum_{\begin{subarray}{c}r\geqslant
0\\\ m\geqslant
1\end{subarray}}^{\infty}\frac{\xi_{r,m}(t_{\lambda}^{(1)})}{n^{\alpha
r+m}}\bigg{\\}}$ $\displaystyle=$
$\displaystyle\omega_{n}^{j}n^{(\alpha+1)\frac{1}{n}}\bigg{\\{}1+\sum_{r=0}^{\lfloor\frac{1}{\alpha}\rfloor}\frac{\xi_{r,1}(\omega_{n}^{j})}{n^{\alpha
r+1}}+O\Big{(}\frac{\log(n)}{n^{2}}\Big{)}\bigg{\\}}$ $\displaystyle=$
$\displaystyle\omega_{n}^{j}\bigg{\\{}1+\sum_{r=0}^{\lfloor\frac{1}{\alpha}\rfloor}\frac{\xi_{r,1}(\omega_{n}^{j})}{n^{\alpha
r+1}}+O\Big{(}\frac{\log(n)}{n^{2}}\Big{)}\bigg{\\}}.$

To reveal the shape in the full expansion we need a third iteration,

$\displaystyle t_{\lambda}^{(3)}$ $\displaystyle\sim$
$\displaystyle\omega_{n}^{j}n^{(\alpha+1)\frac{1}{n}}\bigg{\\{}1+\sum_{\begin{subarray}{c}r\geqslant
0\\\ m\geqslant
1\end{subarray}}^{\infty}\frac{\xi_{r,m}(t_{\lambda}^{(2)})}{n^{\alpha
r+m}}\bigg{\\}}$ $\displaystyle=$
$\displaystyle\omega_{n}^{j}n^{(\alpha+1)\frac{1}{n}}\bigg{\\{}1+\sum_{r=0}^{\lfloor\frac{2}{\alpha}\rfloor}\frac{\xi_{r,1}(\omega_{n}^{j})}{n^{\alpha
r+1}}+\sum_{r=0}^{\lfloor\frac{1}{\alpha}\rfloor}\sum_{\ell=0}^{1}\zeta_{l,r}^{(3)}(\omega_{n}^{j})\frac{\log^{\ell}(n)}{n^{\alpha
r+2}}+O\Big{(}\frac{\log(n)}{n^{3}}\Big{)}\bigg{\\}}$ $\displaystyle=$
$\displaystyle\omega_{n}^{j}\bigg{\\{}1+\sum_{r=0}^{\lfloor\frac{2}{\alpha}\rfloor}\frac{\zeta_{r}^{(4)}(\omega_{n}^{j})}{n^{\alpha
r+1}}+\sum_{r=0}^{\lfloor\frac{1}{\alpha}\rfloor}\sum_{\ell=0}^{1}\zeta_{j,r}^{(5)}(\omega_{n}^{j})\frac{\log^{\ell}(n)}{n^{\alpha
r+2}}+O\Big{(}\frac{\log(n)}{n^{3}}\Big{)}\bigg{\\}},$

where the functions $\zeta_{l,r}^{(3)}$, $\zeta_{r}^{(4)}$, and
$\zeta_{j,r}^{(5)}$ belong to $C^{\infty}(W\setminus\overline{\mathbb{D}})$.
By continuing this process we finish the theorem. ∎

###  2.1 Algorithm

As in previous works, one key observation is that $\omega_{n}^{j}$ can be set
unchanged under several combinations of $n$ and $j$, for example

$\omega_{n}^{j}=\omega_{cn}^{cj},$

for any constant $c$, thus we take $n_{\ell}\equiv\ell n_{1}$ and
$j_{\ell}\equiv\ell j_{1}$ for $\ell\in\mathbb{N}$. The index $j_{\ell}$
depends on $j_{1}$, and similarly, the matrix sizes $n_{\ell}$ depend on
$n_{1}$, but for notation simplicity, we suppressed those dependencies.
Therefore we have

$\sigma_{j_{1}}\equiv\omega_{n_{1}}^{j_{1}}=\omega_{n_{2}}^{j_{2}}=\cdots=\omega_{n_{k-1}}^{j_{k-1}},$

see Figure 2. The sequence $\\{\omega_{n}^{j}\\}_{j=0}^{n-1}$ corresponds to
the $n$ roots of the unity in the complex plane $\mathbb{C}$ and is a regular
partition of the unit circle $\mathbb{T}$ with step angle of $\frac{2\pi}{n}$.
We assume that $n_{1}$ and $k$ are fixed natural numbers and $n\gg n_{1}$.

Our algorithm is composed by three phases: precomputing, extrapolation, and
interpolation. We now describe them in detail.

Figure 2: For $n_{1}=12$ the regular circular grids
$\\{\omega_{n_{\ell}}^{j}\\}$ for $j=1,\ldots,n_{\ell}$ and $\ell=1,2,3$. In
the precomputing phase we will need to calculate the eigenvalues of
$T_{n_{\ell}}(a)$ for $\ell=1,\ldots,k-1$ corresponding to the blue and red
dots combined, but in the interpolation phase, we will use only the
eigenvalues corresponding to the red dots, that is
$\lambda_{j_{\ell}}(T_{n_{\ell}}(a))$ for $j_{1}=1,\ldots,n_{1}$ and
$\ell=1,\ldots,k-1$.

Precomputing. In this phase we need to calculate the eigenvalues of $T_{n}(a)$
for $n=n_{1},\ldots,n_{k-1}$. The numerical computations for the “exact”
eigenvalues are conducted with the BigFloat data type in Julia v.1.7.2 using
1024 bit precision which approximately corresponds to 300 decimal digits of
accuracy. We use an AMD Epyc with 32 cores/64 threads and 1024 GB RAM. The
remaining computations are carried using Mathematica v13.

Extrapolation. Here we need to simplify the expansion in Theorem 2.1, for that
purpose, let $\xi_{1}(n),\xi_{2}(n),\ldots$, be the asymptotic ordering of the
set

$\Big{\\{}\frac{\log^{\ell}(n)}{n^{\alpha r+m}}\colon r\geqslant 0,\
m\geqslant 1,\ \ell=0,\ldots,m-1\Big{\\}}.$

This ordering depends on $\alpha$, for example for $\alpha=\frac{3}{4}$ we
obtain

$\displaystyle\xi_{1}(n)$ $\displaystyle=\frac{1}{n},$
$\displaystyle\xi_{2}(n)$ $\displaystyle=\frac{1}{n^{\alpha+1}},$
$\displaystyle\xi_{3}(n)$ $\displaystyle=\frac{\log(n)}{n^{2}},$
$\displaystyle\xi_{4}(n)$ $\displaystyle=\frac{1}{n^{2}}$
$\displaystyle\xi_{5}(n)$ $\displaystyle=\frac{1}{n^{2\alpha+1}},$
$\displaystyle\xi_{6}(n)$ $\displaystyle=\frac{\log(n)}{n^{\alpha+2}},$
$\displaystyle\xi_{7}(n)$ $\displaystyle=\frac{1}{n^{\alpha+2}},$
$\displaystyle\xi_{8}(n)$ $\displaystyle=\frac{\log^{2}(n)}{n^{3}}.$

For each positive entire number $k$, rewrite the expansion as follows

$t_{j,n}=\omega_{n}^{j}n^{(\alpha+1)\frac{1}{n}}\Big{\\{}1+\sum_{s=1}^{k-1}c_{s}(\omega_{n}^{j})\xi_{s}(n)+E_{k,j,n}\Big{\\}},$
(2.6)

where $E_{k,j,n}=O(\xi_{k}(n))$ is a remainder or error term satisfying
$|E_{k,j,n}|\leqslant\kappa\,\xi_{k}(n)$ for some constant $\kappa$ depending
only on the generating function $a$. For each fixed $j_{1}=1,\ldots,n_{1}$ let
$\sigma_{j_{1}}\equiv\omega_{n_{1}}^{j_{1}}=\cdots=\omega_{n_{k-1}}^{j_{k-1}}$
(see the radially aligned red dots in Figure 2), and apply $k-1$ times the
expansion (2.6), obtaining

$\displaystyle t_{j_{1},n_{1}}$ $\displaystyle=$
$\displaystyle\sigma_{j_{1}}n_{1}^{(\alpha+1)\frac{1}{n_{1}}}\Big{\\{}1+c_{1}(\sigma_{j_{1}})\xi_{1}(n_{1})+\cdots+c_{k-1}(\sigma_{j_{1}})\xi_{k-1}(n_{1})+E_{k,j_{1},n_{1}}\Big{\\}},$
$\displaystyle t_{j_{2},n_{2}}$ $\displaystyle=$
$\displaystyle\sigma_{j_{1}}n_{2}^{(\alpha+1)\frac{1}{n_{2}}}\Big{\\{}1+c_{1}(\sigma_{j_{1}})\xi_{1}(n_{2})+\cdots+c_{k-1}(\sigma_{j_{1}})\xi_{k-1}(n_{2})+E_{k,j_{2},n_{2}}\Big{\\}},$
$\displaystyle\vdots$ $\displaystyle t_{j_{k-1},n_{k-1}}$ $\displaystyle=$
$\displaystyle\sigma_{j_{1}}n_{k-1}^{(\alpha+1)\frac{1}{n_{k-1}}}\Big{\\{}1+c_{1}(\sigma_{j_{1}})\xi_{1}(n_{k-1})+\cdots+c_{k-1}(\sigma_{j_{1}})\xi_{k-1}(n_{k-1})+E_{k,j_{k-1},n_{k-1}}\Big{\\}}.$

Let $\hat{c}_{s}(\sigma_{j_{1}})$ be the approximation of
$c_{s}(\sigma_{j_{1}})$ obtained by removing all the remainder terms
$E_{k,j_{\ell},n_{\ell}}$ and solving the resulting linear system:

$\begin{bmatrix}\xi_{1}(n_{1})&\xi_{2}(n_{1})&\cdots&\xi_{k-1}(n_{1})\\\\[3.44444pt]
\xi_{1}(n_{2})&\xi_{2}(n_{2})&\cdots&\xi_{k-1}(n_{2})\\\\[3.44444pt]
\vdots&\vdots&\ddots&\vdots\\\\[3.44444pt]
\xi_{1}(n_{k-1})&\xi_{2}(n_{k-1})&\cdots&\xi_{k-1}(n_{k-1})\end{bmatrix}\begin{bmatrix}\hat{c}_{1}(\sigma_{j_{1}})\\\\[3.44444pt]
\hat{c}_{2}(\sigma_{j_{1}})\\\\[3.44444pt] \vdots\\\\[3.44444pt]
\hat{c}_{k-1}(\sigma_{j_{1}})\end{bmatrix}=\sigma_{j_{1}}^{-1}\begin{bmatrix}n_{1}^{-(\alpha+1)\frac{1}{n_{1}}}t_{j_{1},n_{1}}\\\
n_{2}^{-(\alpha+1)\frac{1}{n_{2}}}t_{j_{2},n_{2}}\\\ \vdots\\\
n_{k-1}^{-(\alpha+1)\frac{1}{n_{k-1}}}t_{j_{k-1},n_{k-1}}\end{bmatrix}\\!\\!-\\!\\!\begin{bmatrix}1\\\\[3.44444pt]
1\\\\[3.44444pt] \vdots\\\\[3.44444pt] 1\end{bmatrix}.$

Thanks to properties (i)–(vi), the function $a$ is bijective, and hence the
value of each $t_{j_{\ell},n_{\ell}}$ can be calculated as
$a^{-1}(\lambda_{j_{\ell}}(T_{n_{\ell}}(a)))$. In general, the latter
computation can be done numerically with any standard root finder (i.e.
FindRoot in Mathematica or fzero in Matlab). As mentioned in previous works, a
variant of this extrapolation strategy was first suggested by Albrecht
Böttcher in [2, §7] and is analogous to the Richardson extrapolation employed
in the context of Romberg integration [22, §3.4].

Interpolation. The objective of this phase is to estimate
$\hat{c}_{s}(\omega_{n}^{j})$, for any index $j\in\\{1,\ldots,n\\}$, any $n$,
and any level $s=1,\ldots,k-1$. If $\omega_{n}^{j}$ coincides with one of the
points in the grid $\\{\omega_{n_{1}}^{1},\ldots,\omega_{n_{1}}^{n_{1}}\\}$,
then we have the approximations $\hat{c}_{s}(\omega_{n}^{j})$ from the
extrapolation phase for free. In any other case, we will do it by
interpolating the data

$(\omega_{n_{1}}^{1},\hat{c}_{s}(\omega_{n_{1}}^{1})),(\omega_{n_{1}}^{2},\hat{c}_{s}(\omega_{n_{1}}^{2})),\ldots,(\omega_{n_{1}}^{n_{1}},\hat{c}_{s}(\omega_{n_{1}}^{n_{1}})),$

and then evaluating the resulting polynomial at $\omega_{n}^{j}$. This
interpolation can be done in many ways, but to avoid spurious oscillations
explained by the Runge phenomenon [11, p. 78], and following the strategy of
the previous works, we decided to do it considering only the $k-s+8$ nodes in
the grid $\\{\omega_{n_{1}}^{1},\ldots,\omega_{n_{1}}^{n_{1}}\\}$ which are
closest to $\omega_{n}^{j}$. Those nodes can be determined uniquely unless
$\omega_{n}^{j}$ is the mid point of two consecutive nodes in the grid, in
which case we can take any of the two possible choices.

Finally, our eigenvalue approximation with $k$ terms, is given by

$\lambda_{j}^{\textsc{inn}(k)}(T_{n}(a))\equiv
a\Big{(}\omega_{n}^{j}n^{(\alpha+1)\frac{1}{n}}\Big{\\{}1+\sum_{s=1}^{k-1}\hat{c}_{s}(\omega_{n}^{j})\xi_{s}(n)\Big{\\}}\Big{)}.$
(2.7)

######  Remark 2.4.

The approximation (2.7) can be tuned in a number of ways. To guarantee
acceptable overall accuracy, the precomputing phase must be done with a
significant amount of precision digits, let us say 60. In the precomputing and
extrapolation phases, we can manage more levels $k$ than in the interpolation
phase, for example 2 or 3 additional levels will improve the final result.
Finally, in the interpolation phase we can increase the number of interpolated
points.

##  3 Extreme eigenvalues

In this section we develop a numerical algorithm to approximate the extreme
eigenvalues of $T_{n}(a)$. In the paper [4] the authors considered the
generating function

$a(t)=\frac{1}{t}(1-t)^{\alpha}f(t),$

satisfying the following properties:

  1. (i)

The function $f$ belongs to $H^{\infty}$ with $f(0)\neq 0$ and for some
$\varepsilon>0$, $f$ has an analytic continuation to the region

$K_{\varepsilon}\equiv B(1,\varepsilon)\setminus\\{x\in\mathbb{R}\colon
1<x<1+\varepsilon\\},$

and is continuous in
$\hat{K}_{\varepsilon}\equiv\overline{B(1,\varepsilon)}\setminus\\{x\in\mathbb{R}\colon
1<x<1+\varepsilon\\}$. Additionally, $f_{\varphi}(x)\equiv f(1+x\hskip
0.6pt\mathrm{e}^{\hskip 0.2pt\mathrm{i}\varphi})$ belongs to the algebra
$C^{2}[0,\varepsilon]$ for each $\varphi\in(-\pi,\pi]$. Here $B(z_{0},r)$ is
the open disk centered at $z_{0}$ with radius $r$.

  2. (ii)

Let $0<\alpha<1$ be a constant and take the argument
$-\alpha\pi<\arg(1-z)^{\alpha}\leqslant\alpha\pi$ when
$-\pi<\arg(1-z)\leqslant\pi$.

  3. (iii)

The range $\mathcal{R}(a)$ is a Jordan curve in $\mathbb{C}$ with
$\operatorname{wind}_{\lambda}(a)=-1$ for each $\lambda\in\mathcal{D}(a)$.

For the eigenvalues $\lambda_{j}(T_{n}(a))$ located in $W_{0}$, it was proved
that

$\lambda_{j}(T_{n}(a))=\frac{\Lambda_{j}^{\alpha}}{(n+1)^{\alpha}}\\{1+\Delta(\Lambda_{j},n)\\},$
(3.1)

where $\Delta(\Lambda_{j},n)=O(n^{-\alpha})$, as $n\to\infty$, represents the
remainder or error term, which satisfies the inequality
$|\Delta(\Lambda_{j},n)|\leqslant\kappa n^{-\alpha}$ for some constant
$\kappa$ depending only on $a$. The numbers $\Lambda_{j}$ are certain roots of
the function

$F(\Lambda)\equiv\frac{2\pi\hskip 0.2pt\mathrm{i}}{\alpha}\hskip
0.6pt\mathrm{e}^{\Lambda}-\int_{C}\hskip 0.6pt\mathrm{e}^{-\Lambda
u}\beta(u)\operatorname{d\\!}u,$

where $C$ is the radial segment in $\mathbb{C}$ with argument
$\pm\frac{3}{4}\pi$ (the sign depending on the sign of the imaginary part of
$\lambda_{j}(T_{n}(a))$) and $\beta(u)\equiv(u^{\alpha}\hskip
0.6pt\mathrm{e}^{\hskip 0.2pt\mathrm{i}\alpha\pi}-1)^{-1}-(u^{\alpha}\hskip
0.6pt\mathrm{e}^{-\hskip 0.2pt\mathrm{i}\alpha\pi}-1)^{-1}$.

In our paper, instead of the exact calculation of the zeros $\Lambda_{j}$, we
use (3.1) as an inspiration to formulate the following conjecture.

######  Conjecture 3.1.

Let $\varepsilon$ be a small positive number, and $a$ be a generating function
satisfying the previous properties (i)–(iii). There exists an entire number
$k\geqslant 1$ such that the eigenvalues $\lambda_{j}(T_{n}(a))$ with
$j=1,\ldots,\lfloor\varepsilon n\rfloor$ or $j=n-\lfloor\varepsilon
n\rfloor+1,\ldots,n$, satisfy the expansion

$\lambda_{j}(T_{n}(a))=\sum_{\ell=0}^{k-1}\frac{q_{\ell}(j)}{(n+1)^{\alpha+\ell}}+E_{k,j,n},$

where

  * •

$E_{k,j,n}$ is the remainder term which satisfies the inequality
$|E_{k,j,n}|\leqslant cn^{-\alpha-k}$ for some constant $c$ depending only on
$\alpha$;

  * •

the coefficients $q_{\ell}$ $(\ell\geqslant 1)$ are continuous functions from
$\mathbb{T}$ to $\mathbb{C}$ depending only on $a$.

###  3.1 Algorithm

For notation simplicity, take $h=\frac{1}{n+1}$. Since Conjecture 3.1 does not
include a factor of the type $\omega_{j}^{n}$ or $\pi jh$, we cannot handle
the grid trick previously used, hence this algorithm is essentially different
from the one in Subsection 2.1 and from the algorithms employed in previous
works.

Precomputing. This phase is similar to the precomputing phase in Subsection
2.1, but we can work with a different number of levels $k$. However, in order
to avoid a heavy notation, we will use the same letter anyway. In this phase
we calculate the eigenvalues of $T_{n}(a)$ for $n=n_{1},\ldots,n_{k-1}$.

Extrapolation. For some fixed index $j_{0}\leqslant n_{1}$ and each
$j=1,\ldots,j_{0}$ or $j=n-j_{0}+1,\ldots,n$, apply $k$ times the expansion in
Conjecture 3.1, reaching

$\displaystyle\lambda_{j}(T_{n_{1}}(a))$ $\displaystyle=$ $\displaystyle
q_{1}(j)h_{1}^{\alpha}+q_{2}(j)h_{1}^{\alpha+1}+\cdots+q_{k}(j)h_{1}^{\alpha+k-1}+E_{k,j,n_{1}},$
$\displaystyle\lambda_{j}(T_{n_{2}}(a))$ $\displaystyle=$ $\displaystyle
q_{1}(j)h_{2}^{\alpha}+q_{2}(j)h_{2}^{\alpha+1}+\cdots+q_{k}(j)h_{2}^{\alpha+k-1}+E_{k,j,n_{2}},$
$\displaystyle\vdots$ $\displaystyle\lambda_{j}(T_{n_{k}}(a))$
$\displaystyle=$ $\displaystyle
q_{1}(j)h_{k}^{\alpha}+q_{2}(j)h_{k}^{\alpha+1}+\cdots+q_{k}(j)h_{k}^{\alpha+k-1}+E_{k,j,n_{k}},$

where $E_{k,j,n}=O(h^{\alpha+k})$ is a remainder term satisfying the
inequality $|E_{k,j,n}|\leqslant c\,h^{\alpha+k}$ for some constant $c$
depending only on $a$. As in the inner eigenvalue case, let
$\hat{q}_{\ell}(j)$ be the approximation of $q_{\ell}(j)$ obtained by removing
all the remainder terms $E_{k,j,n_{\ell}}$ and solving the resulting linear
system:

$\begin{bmatrix}h_{1}^{\alpha}&h_{1}^{\alpha+1}&\cdots&h_{1}^{\alpha+k-1}\\\
h_{2}^{\alpha}&h_{2}^{\alpha+1}&\cdots&h_{2}^{\alpha+k-1}\\\
\vdots&\vdots&\ddots&\vdots\\\
h_{k}^{\alpha}&h_{k}^{\alpha+1}&\cdots&h_{k}^{\alpha+k-1}\end{bmatrix}\begin{bmatrix}\hat{q}_{1}(j)\\\
\hat{q}_{2}(j)\\\ \vdots\\\
\hat{q}_{k}(j)\end{bmatrix}=\begin{bmatrix}\lambda_{j}(T_{n_{1}}(a))\\\
\lambda_{j}(T_{n_{2}}(a))\\\ \vdots\\\
\lambda_{j}(T_{n_{k}}(a))\end{bmatrix}.$ (3.2)

Finally, our extreme eigenvalue approximation with $k$ terms is given by

$\lambda_{j}^{\textsc{ext}(k)}(T_{n}(a))\equiv\sum_{\ell=0}^{k-1}\frac{\hat{q}_{\ell}(j)}{(n+1)^{\alpha+\ell}}.$
(3.3)

##  4 Numerical experiments

In this section we test the accuracy of the proposed algorithms related to the
formula (2.7) and to the formula (3.3), for approximating the inner and
extreme eigenvalues of $T_{n}(a)$, respectively. Our numerical experiments are
performed by using Mathematica v.12 (64 bit) on a platform with 16GB RAM,
using an Intel processor QuadCore IntelCore i7 2.6 GHz. For a generating
function $a$ and a level $k$, consider the individual absolute errors

$\textrm{AE}_{j,n}^{\textsc{inn}(k)}\equiv|\lambda_{j}(T_{n}(a))-\lambda_{j}^{\textsc{inn}(k)}(T_{n}(a))|,\quad\textrm{AE}_{j,n}^{\textsc{ext}(k)}\equiv|\lambda_{j}(T_{n}(a))-\lambda_{j}^{\textsc{ext}(k)}(T_{n}(a))|.$

To improve the readability of the obtained errors, we decide to define the
maximum error

$\textrm{AE}_{n}^{\textsc{inn}(k)}\equiv\max\\{\textrm{AE}_{j,n}^{\textsc{inn}(k)}\colon
j=\lfloor\varepsilon n\rfloor,\ldots,n-\lfloor\varepsilon n\rfloor\\},$

and the individual relative errors

$\textrm{RE}_{j,n}^{\textsc{ext}(k)}\equiv\frac{\textrm{AE}_{j,n}^{\textsc{ext}(k)}}{|\lambda_{j}(T_{n}(a))|},\qquad\textrm{RE}_{j,n}^{\textsc{inn}(k)}\equiv\frac{\textrm{AE}_{j,n}^{\textsc{inn}(k)}}{|\lambda_{j}(T_{n}(a))|}.$

Given an $\varepsilon>0$, our global individual relative error is

$\textrm{RE}_{j,n}^{\varepsilon}\equiv\begin{cases}\textrm{RE}_{j,n}^{\textsc{ext}(k_{2})}&\quad\mbox{if}\quad
j\in[1,\lfloor\varepsilon n\rfloor]\cup[n-\lfloor\varepsilon n\rfloor+1,n];\\\
\textrm{RE}_{j,n}^{\textsc{inn}(k_{1})}&\quad\mbox{if}\quad
j\in[\lfloor\varepsilon n\rfloor+1,n-\lfloor\varepsilon n\rfloor],\end{cases}$

where $k_{1}$ and $k_{2}$ are the largest available. Note that these constants
are limited by the precomputing phase, and that $\varepsilon$ must be selected
within the range $0<\varepsilon<\frac{n_{1}}{n}$.

For $\alpha\in\mathbb{R}_{+}\setminus\mathbb{Z}$, consider the generating
function

$a(t)=\frac{1}{t}(1-t)^{\alpha}.$

A detailed description of this function and all the related neighborhoods
appears in Figure 1. Using the binomial theorem we obtain

$a(t)=\sum_{j=-1}^{\infty}(-1)^{j+1}\binom{\alpha}{j+1}t^{j},$

then the Fourier coefficients $a_{j}$ are given by
$a_{j}=(-1)^{j+1}\binom{\alpha}{j+1}$ for $j\geqslant-1$ and $a_{j}=0$ in any
other case. It is also worth noticing that since the matrices $T_{n}(a)$ are
all real, then the eigenvalues appear as conjugate pairs and this explains the
symmetry of the error in the related figures.

$n$ | $256$ | $512$ | $1024$ | $2048$ | $4096$ | $8192$
---|---|---|---|---|---|---
$\textrm{AE}_{n}^{\textsc{inn}(1)}$ | $1.5783\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $8.0820\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $4.0956\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $2.0631\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $1.0358\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $5.1906\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$
$\xi_{1}(n)\textrm{AE}_{n}^{\textsc{inn}(1)}$ | $4.0405\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $4.1380\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $4.1938\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $4.2252\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $4.2426\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $4.2521\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$
$\textrm{AE}_{n}^{\textsc{inn}(2)}$ | $3.0955\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $8.8936\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $2.5152\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $7.0464\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $1.9630\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $5.4512\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$
$\xi_{2}(n)\textrm{AE}_{n}^{\textsc{inn}(2)}$ | $5.0716\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $4.9012\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $4.6623\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $4.3934\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $4.1167\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $3.8452\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$
$\textrm{AE}_{n}^{\textsc{inn}(3)}$ | $2.3674\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $6.6709\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $1.8445\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $5.0358\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $1.3625\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $3.6613\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$
$\xi_{3}(n)\textrm{AE}_{n}^{\textsc{inn}(3)}$ | $2.7980\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $2.8032\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $2.7904\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $2.7702\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $2.7482\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $2.7268\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$
$\textrm{AE}_{n}^{\textsc{inn}(4)}$ | $7.8781\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $2.0791\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $5.3823\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $1.3770\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $3.4968\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $8.8421\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$
$\xi_{4}(n)\textrm{AE}_{n}^{\textsc{inn}(4)}$ | $5.1630\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $5.4503\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $5.6438\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $5.7754\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $5.8666\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $5.9338\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$
$\textrm{AE}_{n}^{\textsc{inn}(5)}$ | $1.1053\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $2.0675\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $3.8406\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $7.1041\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $1.3076\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $2.3356\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$
$\xi_{5}(n)\textrm{AE}_{n}^{\textsc{inn}(5)}$ | $1.1589\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $1.2264\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $1.2887\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $1.3484\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $1.4040\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $1.4187\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$
$\textrm{AE}_{n}^{\textsc{inn}(6)}$ | $8.2836\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $1.4196\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $2.3767\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $3.9140\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $6.3845\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.0800\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$
$\xi_{6}(n)\textrm{AE}_{n}^{\textsc{inn}(6)}$ | $6.2656\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $6.4210\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $6.5083\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $6.5549\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $6.5934\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$ | $6.9258\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{0}$
Table 1: The maximum errors $\textup{AE}_{n}^{\textsc{inn}(k)}$, and maximum
normalized errors $\xi_{k}(n)\textup{AE}_{n}^{\textsc{inn}(k)}$ for the levels
$k=1,\ldots,6$, indexes $j\in\big{(}\frac{n}{8},\frac{7n}{8}\big{)}$, and
different matrix sizes $n$, corresponding to the matrix $T_{n}(a)$ with
$a(t)=\frac{1}{t}(1-t)^{\frac{3}{4}}$. We used a grid of size $n_{1}=100$.

Figures 3 and 4 show a comparison between the individual relative errors
$\textrm{RE}_{j,n}^{\textsc{ext}(k)}$ and
$\textrm{RE}_{j,n}^{\textsc{inn}(k)}$, obtained by our expansions for
different number of terms $k$. We can see that our algorithms are producing
smaller errors as $k$ increases, and that the overall reached precision is
very high.

Table 1 shows the maximum individual and normalized inner eigenvalue errors.
We considered only the indexes $j\in\big{(}\frac{n}{8},\frac{7n}{8}\big{)}$,
in this way we are excluding the eigenvalues belonging to the neighborhood
$W_{0}$ corresponding to $\varepsilon=\frac{1}{8}$. The almost constant
behavior of the normalized errors is a verification of the expansion in
Theorem 2.1. We remark that the error $\textup{AE}^{\textsc{inn}(1)}_{n}$ has
order $O(\frac{1}{n})$ which is the same that the distance between two
consecutive eigenvalues have, which means that the approximation with one term
can be only used to understand the distribution of the eigenvalues.

Figure 3: The base-10 logarithm extreme eigenvalue relative errors
$\textup{RE}^{\textsc{ext}(k)}_{j,n}$ $(j=1,\ldots,100)$ for the matrix
$T_{n}(a)$ with $a(t)=\frac{1}{t}(1-t)^{\frac{3}{4}}$, a grid size
$n_{1}=100$, a matrix size $n=8192$, and different number of terms $k=1$
(blue), $k=2$ (red), $k=3$ (black), and $k=6$ (green). All the errors are
plotted over the grid $\big{\\{}\frac{2\pi j}{n+1}\big{\\}}_{j=1}^{n}$. Figure
4: The same as Figure 3 for the inner eigenvalue relative errors
$\textup{RE}_{j,n}^{\textsc{inn}(k)}$ with $\varepsilon=\frac{1}{64}$, and for
$k=3$ (blue), $k=4$ (red), and $k=5$ (green). All the errors are plotted over
the grid $\big{\\{}\frac{2\pi j}{n+1}\big{\\}}_{j=1}^{n}$. Figure 5: The
base-10 logarithm for the global individual relative errors
$\textup{RE}_{j,n}^{\varepsilon}$ for the matrix $T_{n}(a)$ with
$a(t)=\frac{1}{t}(1-t)^{\frac{3}{4}}$, a grid size $n_{1}=100$, levels
$k_{1}=7$, $k_{2}=7$, and different matrix sizes and epsilons: $n=2048$,
$\varepsilon=\frac{1}{30}$ (blue), $n=4096$, $\varepsilon=\frac{1}{45}$ (red),
and $n=8192$, $\varepsilon=\frac{1}{90}$ (green). All the errors are plotted
over the grid $\big{\\{}\frac{2\pi j}{n+1}\big{\\}}_{j=1}^{n}$.

Figure 5 shows the global individual relative errors
$\textrm{RE}_{j,n}^{\varepsilon}$ for different matrix sizes $n$ and epsilons.
We can see the high precision attained by our global algorithm, specially at
the extreme eigenvalues.

##  5 Future work

The case $m<\alpha<m+1$ for $m=1,2,\ldots$ is very interesting for the
applications. We will consider those cases in a future project.

It would be of interest to use such results also for the fast solution of the
resulting large linear systems, so covering the preconditioned case and the
setting of FDEs defined on $d$-dimensional domains with $d\geqslant 2$.

## Acknowledgments

Part of the numerical computations were carried in the computer center Jürgen
Tischer of the mathematics department at Universidad del Valle.

The research of S.M. Grudsky was supported by CONACYT (Mexico) project
“Ciencia de Frontera” FORDECYT-PRONACES/61517/2020 and by Regional
Mathematical Center of the Southern Federal University with the support of the
Ministry of Science and Higher Education of Russia, Agreement
075-02-2021-1386.

The research of the third author is partly supported by the Italian INdAM-GNCS
agency.

## References

  * [1] Bogoya, M., Böttcher, A., and Grudsky, S.M. Asymptotics of individual eigenvalues of a class of large Hessenberg Toeplitz matrices.  Oper. Theory: Adv. Appl. 220 (2012), 77–95.
  * [2] Bogoya, M., Böttcher, A., Grudsky, S.M., and Maximenko, E.A. Eigenvalues of Hermitian Toeplitz matrices with smooth simple-loop symbols.  J. Math. Anal. Appl. 422 (2015), 1308–1334.
  * [3] Bogoya, M., Ekström, S.-E., and Serra-Capizzano, S. Fast Toeplitz eigenvalue computations joining interpolation-extrapolation matrix-less algorithms and simple-loop theory.  Numer. Alg. (2022).
  * [4] Bogoya, M., Grudsky, S.M., and Malysheva, I.S. Extreme individual eigenvalues for a class of large Hessenberg Toeplitz matrices.  Oper. Theory Adv. Appl. 271 (2018), 119–143.
  * [5] Bogoya, M., Grudsky, S.M., Mazza, M., and Serra-Capizzano, S. On the extreme eigenvalues and asymptotic conditioning of a class of Toeplitz matrix-sequences arising from fractional problems.  Linear Multilin. Alg. (2022).
  * [6] Bogoya, M., Grudsky, S.M., Serra-Capizzano, S., and Tablino-Possio, C. Fine spectral estimates with applications to the optimally fast solution of large FDE linear systems.  BIT Numer. Math. (2022).
  * [7] Bogoya, M., and Serra-Capizzano, S. Eigenvalue superposition expansion for Toeplitz matrix-sequences, generated by linear combinations of matrix-order dependent symbols, and applications to fast eigenvalue computations.  arXiv:2112.11794 (2022).
  * [8] Bogoya, M., Serra-Capizzano, S., and Vassalos, P. Fast Toeplitz eigenvalue computations joining interpolation-extrapolation matrix-less algorithms and simple-loop theory: The preconditioned setting.  arXiv:2203.11338 (2022).
  * [9] Böttcher, A., and Grudsky, S.M. On the condition numbers of large semi-definite Toeplitz matrices.  Linear Algebra Appl. 279, 1/3 (1998), 285–301.
  * [10] Böttcher, A., and Grudsky, S.M. Spectral properties of banded Toeplitz matrices.  Society for Industrial and Applied Mathematics (SIAM), 2005.
  * [11] Davis, P.J. Interpolation and approximation.  Dover, New York, 1975.
  * [12] Donatelli, M., Garoni, C., Manni, C., Serra-Capizzano, S., and Speleers, H. Robust and optimal multi-iterative techniques for IgA Galerkin linear systems.  Comput. Methods Appl. Mech. Engrg. 284 (2015), 230–264.
  * [13] Donatelli, M., Mazza, M., and Serra-Capizzano, S. Spectral analysis and structure preserving preconditioners for fractional diffusion equations.  J. Comput. Phys. 307 (2016), 262–279.
  * [14] Ekström, S.-E., Furci, I., and Serra-Capizzano, S. Exact formulae and matrix-less eigensolvers for block banded Toeplitz-like matrices.  BIT 58, 4 (2018), 937–968.
  * [15] Ekström, S.-E., and Garoni, C. A matrix-less and parallel interpolation-extrapolation algorithm for computing the eigenvalues of preconditioned banded symmetric Toeplitz matrices.  Numer. Alg. 80 (2019), 819–848.
  * [16] Ekström, S.-E., Garoni, C., and Serra-Capizzano, S. Are the eigenvalues of banded symmetric Toeplitz matrices known in almost closed form?  Exper. Math. 27, 4 (2018), 478–487.
  * [17] Fedoryuk, M.V. The saddle-point method [Russian].  Nauka, Moscow, 1977.
  * [18] Parter, S.V. Extreme eigenvalues of Toeplitz forms and applications to elliptic difference equations.  Trans. Amer. Math. Soc. 99 (1961), 153–192.
  * [19] Parter, S.V. On the extreme eigenvalues of Toeplitz matrices.  Trans. Amer. Math. Soc. 100 (1961), 263–276.
  * [20] Serra-Capizzano, S. On the extreme eigenvalues of Hermitian (block) Toeplitz matrices.  Linear Algebra Appl. 270 (1998), 109–129.
  * [21] Serra-Capizzano, S., and Tilli, P. Extreme singular values and eigenvalues of non-Hermitian block Toeplitz matrices.  J. Comput. Appl. Math. 108, 1/2 (1999), 113–130.
  * [22] Stoer, J., and Bulirsch, R. Introduction to numerical analysis, third ed.  Springer, 2010.
  * [23] Tilli, P. Universal bounds on the convergence rate of extreme Toeplitz eigenvalues.  Linear Algebra Appl. 366 (2003), 403–416.
  * [24] Widom, H. On the eigenvalues of certain Hermitian operators.  Trans. Amer. Math. Soc. 88 (1958), 491–522.