Eigenvalue Clusters of Large Tetradiagonal Toeplitz Matrices

Toeplitz matrices are typically non-Hermitian and hence they evade the well-elaborated machinery one can employ in the Hermitian case. In a pioneering paper of 1960, Palle Schmidt and Frank Spitzer showed that the eigenvalues of large banded Toeplitz matrices cluster along a certain limiting set which is the union of finitely many closed analytic arcs. Finding this limiting set nevertheless remains a challenge. We here present an algorithm in the spirit of Richard Beam and Robert Warming that reduces testing O(N2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(N^2)$$\end{document} points in the plane for membership in the limiting set by testing only O(N) points along a one-dimensional curve. For tetradiagonal Toeplitz matrices, we describe all types of the limiting sets, we classify their exceptional points, and we establish asymptotic formulas for the analytic arcs near their endpoints.


Introduction
It is well known that the properties of truncated n × n Toeplitz matrices as n goes to infinity are captured by the series a(z) = ∞ j=−∞ a j z j (z ∈ C). We therefore denote the matrix (1.1) by T n (a). The first Szegő limit theorem [12] describes the asymptotic eigenvalue distribution of T n (a) in the case where the restriction of a to the complex unit circle T is (the Fourier series of) a real-valued L ∞ function. More recent studies extend this to realvalued L 1 functions and even more general Hermitian situations; see, e.g., [6,[14][15][16]. The problem is of a completely different nature if the restriction a|T is complex-valued. Here the path-breaking result is due to the 1960 paper [10] by Schmidt  Throughout what follows we always assume that a j ∈ C, r ≥ 1, s ≥ 1, a −r = 0, and a s = 0. Schmidt and Spitzer [10] showed that the set σ(T n (a)) of the eigenvalues of T n (a) converges in the Hausdorff metric to some limiting set Λ(a) which is the union of finitely many analytic arcs and their endpoints and which does not contain isolated points. Ullman [17] later proved that Λ(a) is connected. It is tempting to say that once we know that the limiting set exists, we could simply compute the eigenvalues of T n (a) for n = 100 and this should give a good approximation to Λ(a). Well, let a(z) = z 3 − 3(1 + i)z 2 + 7iz + 4(1 − i) − 2z −1 . (1.2) We thus have a pentadiagonal matrix. Figure 1a shows the curve a(T) and the eigenvalues of T 100 (a) delivered by Matlab. Clearly, part of the eigenvalues are erroneous. Even after zooming in it is impossible to guess what Λ(a) is. Already Schmidt and Spitzer employed the key trick with banded Toeplitz matrices: namely, if D = diag( j−1 ) n j=1 , then D (a j−k ) n j,k=1 D −1 = ( j−k a j−k ) n j,k=1 (1.3) and hence T n (a) = (a j−k ) n j,k=1 and ( j−k a j−k ) n j,k=1 have the same eigenvalues. The last matrix may be written as T n (a ) with a (z) = s j=−r j a j z j . Figure 1b-d shows a (T) and the (in part erroneous) eigenvalues of T 100 (a ) computed with Matlab for three choices of the parameter . Without knowing more on Λ(a) it is difficult to make a reliable guess.
Let us turn to the tetradiagonal case. A subset of C consisting of two distinct points μ, ν ∈ C and an analytic arc (without self-intersection) joining these two points will be denoted by [μ ∼ ν]. Part of our main result is as follows.
Let γ be the piece of the curve {b ∈ C : |1 + b| = 2|b| 2 } that lies in the closed disk {b ∈ C : |1 + b| ≤ 1} (see Fig. 7) and let Γ be the curve  The ± indicates that we take both values of (1 + b + b 2 ) 3/2 . The curve Γ is the boundary of the blue domain Ω shown in Fig. 2. The two red points are ±3 √ 3 and they correspond to b = −1/2 ∈ γ. As will be seen below, the key trick allows us to assume without loss of generality that a(z) = z 2 + cz + cz −1 with c ∈ C\{0}. Denote the roots of a (z) = 0 by t 1 , t 2 , t 3 and order them so that |t 3 | ≤ |t 2 | ≤ |t 1 |. Put λ k = a(t k ). Then Λ(a) is a set of one of the following types. If c is in the interior of Ω, then the three points λ 1 , λ 2 , λ 3 are distinct and   Schmidt and Spitzer [10] gave two descriptions of the limiting set Λ(a). First, they showed that where σ(T (a )) is the spectrum of the Toeplitz operator generated by the infinite Toeplitz matrix ( j−k a j−k ) ∞ j,k=1 on 2 . Note that σ(T (a )) is just the union of the curve a (T) and the points in the plane about which this curve has nonzero winding number. For the second description, pick λ ∈ C, consider the equation a(z) = λ (which after multiplication by z r becomes an equation with an algebraic polynomial of degree r + s), and order the solutions z j (λ) (j = 1, . . . , r + s) of this equation so that The original reasoning of Schmidt and Spitzer was significantly simplified and accompanied with the density function for the eigenvalues along Λ(a) by Hirschman [7] and Widom [18,19]; see also Chapter 11 of [2]. A new view at the matter is given by Duits and Kuijlaars in [5]. We should also note that in the tridiagonal case a(z) = a 1 z + a 0 + a −1 z −1 we have If a(z) = Recently, Shapiro andŠtampach [11] very beautifully characterized the Laurent polynomials for which Λ(a) is a subset of the real line R: this happens if and only if the pre-image a −1 (R) contains a Jordan curve. We refer to the books [2,3,6,8] for more on the eigenvalue cluster problem for Toeplitz matrices. An analogue of the Schmidt-Spitzer result for Wiener-Hopf integral operators is in [4]. We also want to mention the message of [9,13] according to which in certain problems the eigenvalues of non-normal operators may drastically mislead us. Finally, recall that T n (a) and T n (a ) have the same eigenvalues. It follows in particular that If a 1 = 0, this reads and taking so that 2 a 2 = −1 a −1 , we obtain that σ(T n (a)) = σ T n (a 2 z 2 + a 0 + a −1 z −1 ) = a 0 + 2 a 2 σ T n (z 2 + z −1 ) , which leads to a concrete single case. If a 1 = 0, we choose so that a 1 = −1 a −1 and get . The limiting set Λ(z 2 + z −1 ) was already determined by Schmidt and Spitzer: it is the star-shaped set Thus, we are left with the one-parameter family T n (z 2 +cz+cz −1 ), c ∈ C\{0}.
If |z r+1,k (ϕ)| = |z r,k (ϕ)|, then λ k (ϕ) ∈ Λ(a), and otherwise λ k (ϕ) / ∈ Λ(a). To tackle the remaining case ϕ = 0, we replace Eq. (2.1) by lim ϕ→0 (a(z) − a(e iϕ z))/ϕ = 0, that is, by a (z) = 0. Let t k (k = 1, . . . , r +s) be the solutions of this equation and put λ k (0) = a(t k ) for each k. Then solve a(z) = λ k (0) for each k and sort the roots z j,k (0) so that In summary, instead of finding the zeros of a polynomial of degree r + s parametrized by λ on a grid of order O(N 2 ) covering a bounded subset of C (one may take the convex hull of the range a(T) or even better some σ(T (a )) ), we have to find the zeros of r + s + 1 polynomials of degree r + s parametrized by ϕ ∈ (−π, π] and thus by points of a grid of order O(N ).
If u satisfies a(u) − a(e iϕ u) = 0, then v = e iϕ u satisfies the equation a(v) − a(e −iϕ v) = 0. Thus, the solutions u 1 (ϕ), . . . , u r+s (ϕ) coincide with the solutions u 1 (−ϕ), . . . , u r+s (−ϕ), which implies that after appropriate sorting we get λ k (ϕ) = λ k (−ϕ) for all k. We therefore may restrict the algorithm to ϕ ∈ [0, π]. Example 2.1. Figure 4 illustrates the algorithm for the septadiagonal matrix T n (a) induced by Example 2.2. We take the pentadiagonal matrix generated by the polynomial (1.2). Recall that the numerically obtained eigenvalues are in Fig. 1. The algorithm applied to a with = 0.4 yields Fig. 5. In the left picture, the numerically computed eigenvalues of T 100 (a ) are shown as blue circles and the limiting set computed with the algorithm is in green. The right picture is a zoom in of the left.
The rest of the paper is devoted to tetradiagonal Toeplitz matrices. In that case we accompany the topological classification of the point in Λ(a) into regular end exceptional points by the following analytical classification. In accordance with the definition given above, we call λ ∈ Λ(a) a simple point if We refer to a point λ ∈ Λ(a) as a branch point if two (or all three) of the points z 1 (λ), z 2 (λ), z 3 (λ) coincide. In that case we may label the points so that The remaining points λ ∈ Λ(a) are the points for which z 1 (λ), z 2 (λ), z 3 (λ) are three distinct points satisfying |z 1 (λ)| = |z 2 (λ)| = |z 3 (λ)|.
These points will be called multiple points.

Branch Points
Let for the moment a(z) = s j=−r a j z j with r, s ≥ 1 and a −r a s = 0 be a general Laurent polynomial. A point λ 0 ∈ Λ(a) is called a branch point if the equation λ 0 = a(z) has a root of multiplicity at least 2. This is equivalent to saying that λ 0 = a(z 0 ) and that z 0 is a multiple zero of the polynomial z r (a(z) − λ 0 ). The latter happens if and only if z 0 is a root of We now return to the tetradiagonal case a(z) = z 2 + cz + cz −1 with c = 0. Then a (z) = 2z 3 + cz 2 − c z 2 has three (possibly coinciding) zeros. We denote these zeros by t 1 , t 2 , t 3 and label them so that |t 1 | ≥ |t 2 | ≥ |t 3 |.
The blue set Ω in Fig. 2 is closed. Its left and right points are ±3 √ 3, and the upper and lower points are ±i. The set Ω is symmetric about the real and the imaginary axes. Thus, it results from one of its quarters lying in one of the four quadrants by reflexions on the real and imaginary axes. Recall the connection between b and c given by Lemma 4.1. In dependence of the choice of σw 3 , the closure of D + is mapped by onto either the quarter in the first quadrant or the quarter in the third quadrant. To be specific, suppose this quarter is the one in the first quadrant. Then the segment from b = −1/2 to A is mapped to the segment from the (red) point 3 √ 3 to the origin, the arc from A to G + is mapped to the segment from 0 to i, and the arc of γ from G + to b = −1/2 is mapped to the boundary piece of Ω from i to 3 √ 3. The other choice −σw 3 maps the closure of D + to the quarter of Ω in the third quadrant. The closure of D − is mapped into the union of the quarters in the second and fourth quadrants.
Combining the description of Ω with Lemma 4.1 we arrive at the following result. Remark 4.4. The parameter c is immediately available but checking whether c is in Ω or on ∂ Ω may be difficult if c is close to ∂ Ω. However, one may proceed as follows. First solve the cubic equation a (z) = 0, denote the roots by t 1 , t 2 , t 3 , label them so that |t 1 | ≥ |t 2 | ≥ |t 3 |, and put b = t 3 /t 1 . Then We now turn to the following question: how many analytic arcs do come out of the branch points and what is their asymptotic behavior near the branch points? Here is the result. these arcs have the following asymptotic representations in a neighborhood of the points λ 2 , λ 3 : (b) If c ∈ int Ω, then there are three analytic arcs 1 , 2 , and 3 belonging to Λ(a) and coming from the points λ 1 = a(t 1 ), λ 2 = a(t 2 ) and λ 3 = a(t 3 ), respectively, and these arcs have an asymptotic representation of the form (4.4) for j = 1, 2, 3. (c) If c ∈ ∂Ω\{±3 √ 3}, then there are two analytic arcs that come from λ 2 = a(t 2 ) and λ 3 = a(t 3 ) with asymptotic representations of the form (4.4) in a neighborhood of these points. The roots of the equation a(z) = λ 1 are t 1 , t 1 , e iϕ0 t 1 and there are two analytic arcs ± 1 that come from λ 1 = a(t 1 ) and have the representations where ϕ ∈ (ϕ 0 , ϕ 0 + ε) for + 1 and ϕ ∈ (ϕ 0 − ε, ϕ 0 ) for − 1 . These two arcs make a cusp at the point λ 1 . Proof. Let λ 0 ∈ Λ(a) be a branch point and let t 0 ∈ {t 1 , t 2 , t 3 } be a point such that a(t 0 ) = λ 0 and a (t 0 ) = 0. The roots of the equation a(z) = λ 0 are t 0 , t 0 , z 3 . From the proof of Lemma 4.2 we know that |t 0 | < |z 3 | if and only if |t 0 | 3 < 2|t 1 t 2 t 3 |. The latter inequality is always true for t 0 ∈ {t 2 , t 3 } and it is satisfied for t 0 = t 1 if and only if |1 + b| < 2|b| 2 , which happens if and only if c ∈ int Ω. We know from Sect. 2 that if λ ∈ Λ(a) is sufficiently close to λ 0 , then λ = a(u(ϕ)) with a function u = u(ϕ) that is implicitly given by the equation a(u) − a(e iϕ u) = 0 and u(0) = t 0 . The solutions of the equation a(z) = λ are u(ϕ), e iϕ u(ϕ), v(ϕ) with u(ϕ) close to t 0 and v(ϕ) close to z 3 . Thus, if |t 0 | < |z 3 |, then any u(ϕ) satisfying a(u) − a(e iϕ u) = 0 and u(0) = t 0 yields an arc starting at t 0 . We may suppose that ϕ ranges over (0, ε) with sufficiently small ε > 0, because, as noted in Sect. 2, a(u(−ϕ)) = a(u(ϕ)).
We look for u(ϕ) in the form Then Taking into account that The coefficients B 1 and B 2 must be zero. If a (t 0 ) = 0, then c = −t 3 0 , and as also a (t 0 ) = 0, we conclude that t 0 = ± √ 3 and hence c = ±3 √ 3. This case will be treated separately in part (d). Thus, let a (t 0 ) = 0. Then the equation B 1 = 0 yields It follows that B 2 equals Replacing (4.7) by u(ϕ) = t 0 + ∞ k=1 d k ϕ k , taking the full power series for e iϕ and a(z), and proceeding as above, one will see that for k ≥ 2 the coefficient d k enters B k linearly and is completely determined by d j for 1 ≤ j ≤ k − 1. Consequently, there is a unique arc starting at λ 0 . Inserting (4.7), (4.9), (4.10) into the expansion we get the expression on the right of (4.4) for the arc. At this point we have proved all assertions concerning λ 2 and λ 3 (because for them |t 0 | < |z 3 |) and we have also proved part (b), since, as said, for c ∈ int Ω we have |t 0 | < |z 3 |, too. We turn to part (c). In that case there is a positive number ϕ 0 such that {t 0 , t 0 , e iϕ0 t 0 } is the solution set of the equation a(z) = λ 0 . We have to look for solutions u = u(ϕ) of the equation a(u) − a(e iϕ u) = 0 satisfying u(0) = t 0 or u(ϕ 0 ) = e iϕ0 t 0 . As noted above, for the initial condition u(0) = t 0 we need only take ϕ from [0, ε), but the initial condition u(ϕ 0 ) = e iϕ0 t 0 requires taking ϕ from (ϕ 0 − ε, ϕ 0 + ε). We first consider the initial condition u(0) = t 0 . As above, we get a unique solution u(ϕ) = t 0 + ∞ k=1 d k ϕ k with d 1 , d 2 given by (4.9), (4.10). Vieta for the equation a (t 0 ) = 6 e iϕ0 t 0 and results in Consequently, It follows that |u(ϕ)| > |t 0 |. The zeros of z(a(z) − λ) = z 3 + cz 2 − λz + c are u(ϕ), e iϕ u(ϕ), v(ϕ), and Vieta tells us that e iϕ u 2 (ϕ)v(ϕ) = −c. Thus, and because λ ∈ Λ(a) would require that |u(ϕ)| = |e iϕ u(ϕ)| ≤ |v(ϕ)| we conclude that (4.7) does not give an arc of Λ(a).
We are left with the equation a(u) − a(e iϕ u) = 0 in neighborhood of the point ϕ = ϕ 0 . Using power series one can show that the solution is unique. We may write and calculations as before give It follows that We emphasize that in the previous lemma we admit ϕ 1 ≤ ϕ 2 ≤ ϕ 3 . The cases ϕ 1 = ϕ 2 and ϕ 2 = ϕ 3 occur for branch points. For multiple points we have strict inequalities.
Let D 2 be the set of all d ∈ C for which Eq. (5.1) has a solution p ∈ [−2, 2] and ψ ∈ R. Note that the function h p : R → C is 2π-periodic and that h −p (R) = h p (R). Figure 8a shows the ranges h p (R) for several values of p and suggests that D 2 is the set bounded by h 2 (R). This can indeed be proved, but we do not need this for our purposes. Moreover, writing −c = r 3 e iθ and d = r 2 e 2πi/3 as in Lemma 5.1, we get The map c → (−c) 2/3 maps Ω to the set shown in Fig. 8b. The blue set is the image of the part of Ω that lies on the right of the imaginary axis for some k ∈ {0, 1, 2}, the red and green sets are the images for the other two values of k. The left piece of Ω is mapped onto the same set. Thus, the set in Fig. 8b is covered twice. Comparison of Fig. 8a, b indicates that D 2 = (−Ω) 2/3 , and since Ω = −Ω, that actually D 2 = Ω 2/3 . This can again be rigorously verified, but we will not embark on this issue here. The conclusion is that Λ(a) has a point λ 0 for which the equation a(z) = λ 0 has three solutions of the same modulus if and only if (−c) 2/3 ∈ D 2 and that this happens if and only if c ∈ Ω. The result of genuine interest to us is that Λ(a) has a multiple point, that is, a point λ 0 such that a(z) = λ 0 is satisfied by three pairwise distinct numbers of the same modulus, if and only if c ∈ int Ω. This will be a consequence of Theorem 6.1 below.
Theorem 5.2. Let λ 0 be a multiple point of Λ(a) and let re iϕ1 , re iϕ2 , re iϕ3 with −π < ϕ 1 < ϕ 2 < ϕ 3 ≤ π be the solutions of the equation a(z) = λ 0 . Then there are exactly three analytic arcs of Λ(a) coming from λ 0 . Their local representations are In (5.8) one has to take ϕ > ϕ 0 or ϕ < ϕ 0 to enforce that Proof. We proceed as in the case of branch points. Now the equation a(z) − λ 0 = 0 has the solutions re iϕ1 , re iϕ2 , re iϕ3 , where r > 0, ϕ j ∈ (−π, π], and ϕ 1 < ϕ 2 < ϕ 3 . In order to find the analytic arcs that pass through λ 0 in a neighborhood around it, we consider again the equation a(u) − a(e iϕ u) = 0.
We have a(u 0 ) − a(e iϕ0 u 0 ) = 0 in the following three cases: Pick one of these cases. In a neighborhood of (u 0 , ϕ 0 ) we look for u(ϕ) in the form Of course, q depends on j and will eventually become q j . We obtain and thus implying that a (u 0 )q − e iϕ0 a (e iϕ0 u 0 )(q + iu 0 ) = 0. This gives (5.10) Let us calculate a (u 0 ) at u 0 = re iϕ1 . We have Using the Vieta equalities (5.3) we get Inserting the two expressions obtained into (5.10) we arrive at We have Re (q/u 0 ) = 0. Indeed, if Re (q/u 0 ) = 0, then, by (5.10), s(u 0 , ϕ 0 ) must be a real number, which implies that Im e i(ϕ2−ϕ1)/2 = 0, whence ϕ 2 = ϕ 1 , contradicting our requirement ϕ 2 > ϕ 1 . From (5.9) we now get It follows that a(u(ϕ)) belongs to Λ(a) if and only if since if we have (5.11), then the roots of the equation a(z) = λ(ϕ) with λ(ϕ) = a(u(ϕ)) have absolute value smaller than r, which means that the third root v(ϕ) of this equation has absolute value |v(ϕ)| = |c| |u(ϕ)| 2 = r 3 |u(ϕ)| 2 > r, while if (5.11) does not hold, then |v(ϕ)| < |u(ϕ)|. The reasoning is analogous for j = 2 and j = 3. Thus, we have shown that three analytic arcs with the asserted local representation start at λ 0 . Inserting the full power series for u(ϕ) in the equation a(u) − a(e iϕ u) = 0 one can check that all its coefficients are uniquely determined, which implies that no more than three analytic arcs come from λ 0 .

The Types of the Limiting Set
We say that Λ(a) is of type L j (j = 1, 2, 3) if Λ(a) is the closure of the union of exactly j analytic arcs and Λ(a) has exactly j + 1 exceptional points. Theorem 6.1. Let a(z) = z 2 + cz + cz −1 with c = 0. Denote the roots of a (z) = 0 by t 1 , t 2 , t 3 and order them so that |t 3 Proof. If c = ±3 √ 3, then Λ(a) is a line segment by Theorem 4.5(d). Let c / ∈ Ω. Then, by Theorems 4.3 and 4.5, Λ(a) has exactly λ 2 and λ 3 as branch points and in each of them there is one outgoing arc. If Λ(a) had a multiple point λ 0 , then due to Theorem 5.2 there would be three arcs starting at λ 0 . Note that, by Lemma 5.1, Λ(a) can have at most one multiple point. It is impossible to weld these five arcs together without creating additional exceptional points and without increasing the number of outgoing arcs in λ 2 and λ 3 . Thus, Λ(a) cannot posses a multiple point and we must have Λ(a) = [λ 2 ∼ λ 3 ]. Now suppose c ∈ int Ω. Then, again by Theorems 4.3 and 4.5, λ 1 , λ 2 , λ 3 are branch points with a single arc starting in each of them. One cannot glue together these three arcs without creating additional exceptional points or increasing the number of outgoing arcs in the branch points. Consequently, Λ(a) must have a multiple point λ 0 with three outgoing arcs. The only possibility to join these arcs in order to get a connected set with no more exceptional points and without increasing the number of outgoing arcs in the branch points is Λ(a) = [λ 0 ∼ λ 1 ] ∪ [λ 0 ∼ λ 2 ] ∪ [λ 0 ∼ λ 3 ]. By Lemma 5.1, λ 0 = −|c| 4/3 .
Finally, let c = ±3 √ 3 belong to ∂ Ω. Theorems 4.3 and 4.5 tell us that λ 1 , λ 2 , λ 3 are branch points and that λ 1 has two outgoing arcs and λ 2 , λ 3 have one outgoing arc. If Λ(a) had a multiple point, it would be impossible to join the seven arcs without getting new exceptional points or increasing the number of outgoing arcs. It follows that Λ(a) has no multiple point and that Λ(a) = [λ 1 ∼ λ 2 ] ∪ [λ 1 ∼ λ 3 ].
The two arcs (4.5) emerging for c ∈ ∂Ω\{±3 √ 3} make a cusp at the point λ 1 . If c moves along ∂Ω to ±3 √ 3, one of these arcs becomes smaller and smaller, and at c = ±3 √ 3 it disappears while the other arc becomes the line segment [−9, 11.25]. This line segment changes to an analytic arc outside Ω whereas inside Ω it transforms into a star with three analytic legs by giving birth to two legs at one of its ends. At c = ±i, the set Λ(a) is formed by two arcs making a cusp. As c moves from ±i into the exterior of Ω, the cusp is bulging out to an analytic curve, and as c moves from ±i into the interior of Ω, the cusp mutates into a star with three legs by creating a new arc at the cusp point and opening the zero angle of the two arcs of the cusp to a positive angle.
We excluded the parameter c = 0. Trivially, in that case Λ(z 2 ) = {0}. If c = 0 approaches the origin, then the star Λ(z 2 + cz + cz −1 ) with the three legs becomes smaller and smaller and eventually it collapses to the point {0}. Proof. Suppose first that ±c ∈ (3 √ 3, ∞). Such c may be represented in the form (4.2) with −1/2 < b < 0. It follows that w = (1 + b + b 2 ) 1/2 may be taken as a positive real number and (4.2) then tells us that all t k are also real numbers. Consequently, λ k = a(t k ) = t 2 k + ct k + ct −1 k ∈ R. The set Λ(a) is the intersection of the spectra of the operators T (a ) for ∈ (0, ∞). Since a (z) = a (z), the curves a (T) are symmetric about the real line. Hence so also are the spectra of all T (a ) and thus also Λ(a). But