## Spectral Approximation for Quasiperiodic Jacobi Operators

Charles Puelz, Mark Embree and Jake Fillman

Abstract. Quasiperiodic Jacobi operators arise as mathematical models of quasicrystals and in more general studies of structures exhibiting aperiodic order. The spectra of these self-adjoint operators can be quite exotic, such as Cantor sets, and their fine properties yield insight into the associated quantum dynamics, that is, the one-parameter unitary group that solves the time-dependent Schr¨dinger equation. Quasiperio odic operators can be approximated by periodic ones, the spectra of which can be computed via two finite dimensional eigenvalue problems. Since long periods are necessary for detailed approximations, both computational efficiency and numerical accuracy become a concern. We describe a simple method for numerically computing the spectrum of a periodK Jacobi operator in O K ( 2 ) operations, then use the algorithm to investigate the spectra of Schr¨dinger operators with Fibonacci, peo riod doubling, and Thue-Morse potentials.

Mathematics Subject Classification (2010). Primary 47B36, 65F15, 81Q10; Secondary 15A18, 47A75.

Keywords. Jacobi operator, Schr¨dinger operator, quasicrystal, Fibonacci, o period doubling, Thue-Morse.

## 1. Introduction

For given sets of parameters { a n } { , b n } ∈ glyph[lscript] ∞ ( ❩ ), the corresponding Jacobi operator J : glyph[lscript] 2 ( ❩ ) → glyph[lscript] 2 ( ❩ ) is defined entrywise by

$$( \mathcal { O } \psi ) _ { n } = a _ { n - 1 } \psi _ { n - 1 } + b _ { n } \psi _ { n } + a _ { n } \psi _ { n + 1 } \quad \quad ( 1. 1 )$$

for all { ψ n } ∈ glyph[lscript] 2 ( ❩ ). When there exists a positive integer K such that

$$a _ { n } = a _ { n + K }, \ \ b _ { n } = b _ { n + K } \quad \text{for all } n \in \mathbb { Z },$$

the Jacobi operator is said to be periodic with period K . Here we are interested in computing the spectrum σ ( J ) of such an operator when the period K is very long, as a route to high fidelity numerical approximations of the

more intricate spectra of operators with aperiodic coefficients. These spectra are important, as they can help us understand the quantum dynamics of solutions to the time-dependent Schr¨dinger equation [31]. o

The spectrum of a periodK Jacobi operator can be calculated from classical Floquet-Bloch theory, the relevant highlights of which we briefly recapitulate for later reference. (Our presentation most closely follows Toda [48, Ch. 4]; see also, e.g., [44, Ch. 5], [46, Ch. 7],[49].) For a scalar E , any solution to J ψ = Eψ satisfies

$$\begin{bmatrix} \psi _ { p K + 1 } \\ \psi _ { p K } \end{bmatrix} = M _ { K } ^ { p } \begin{bmatrix} \psi _ { 1 } \\ \psi _ { 0 } \end{bmatrix}$$

for each p ∈ ❩ , where M K ≡ M K ( E ) denotes the 2 × 2 monodromy matrix

$$M _ { K } = \begin{bmatrix} \frac { E - b _ { K } } { a _ { K } } & - \frac { a _ { K - 1 } } { a _ { K } } \\ 1 & 0 \end{bmatrix} \cdots \begin{bmatrix} \frac { E - b _ { 2 } } { a _ { 2 } } & - \frac { a _ { 1 } } { a _ { 2 } } \\ 1 & 0 \end{bmatrix} \begin{bmatrix} \frac { E - b _ { 1 } } { a _ { 1 } } & - \frac { a _ { K } } { a _ { 1 } } \\ 1 & 0 \end{bmatrix}. \quad \text{ (1.3)}$$

Note that

$$\det _ { \substack { M _ { K } \\ \dots \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \dots \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \                                                                                                                                                                                                        
  \det _ { \substack { M _ { K } \\ \dots \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, -. \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \. \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \: - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, -. \, - \, - \, - \, - \, -. \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - - \, - \, - \, - \, - \, -. \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \, - \sigma _ { 1 } }$$

Now when J is periodic, E ∈ σ ( J ) provided J ψ = Eψ with bounded nontrivial ψ = { ψ n } , 1 which by (1.2) and (1.4) requires the eigenvalues of M K to have unit modulus. From

$$\det ( \gamma - M _ { K } ) = \gamma ^ { 2 } - \text{tr} ( M _ { K } ) \gamma + 1,$$

where tr( ) denotes the trace, the eigenvalues of · M K have unit modulus when

$$- \, 2 \leq \text{tr} ( M _ { K } ) \leq 2.$$

Since tr( M k ) is a degreeK polynomial in E , in principle we can find σ ( J ) by solving the two polynomial equations tr( M K ( E )) = ± 2, giving σ ( J ) as the union of K real intervals. For large K such numerical computations can incur significant errors, as illustrated in [14, Sec. 7.1]. Alternatively, note that if [ ψ , ψ 1 0 ] T is an eigenvector of M K associated with eigenvalue γ = e iθ , then by definition (1.1) and periodicity, J ψ = Eψ and ψ K j + = γψ j imply

$$\text{by comiton} \, ( 1. 1 ) \, \text{and peroicity}, \, \partial \psi = \mathcal { L } \psi \, \text{and} \, \psi _ { K + j } = \gamma \psi _ { j } \, \text{imply} \\ \begin{bmatrix} b _ { 1 } & a _ { 1 } & & e ^ { - i \theta } a _ { K } \\ a _ { 1 } & b _ { 2 } & \ddots & \\ & \ddots & \ddots & \ddots \\ e ^ { i \theta } a _ { K } & & a _ { K - 1 } & a _ { K - 1 } & b _ { K } \end{bmatrix} \begin{bmatrix} \psi _ { 1 } \\ \psi _ { 2 } \\ \vdots \\ \psi _ { K - 1 } \\ \psi _ { K } \end{bmatrix} = \mathcal { E } \begin{bmatrix} \psi _ { 1 } \\ \psi _ { 2 } \\ \vdots \\ \psi _ { K - 1 } \\ \psi _ { K } \end{bmatrix}, \\ \text{where all unspecified entries in the matrix equal zero. Solving} \, \text{$this $K\times K$}$$

where all unspecified entries in the matrix equal zero. Solving this K × K symmetric matrix eigenvalue problem for any θ ∈ [0 , 2 π ) gives K points in σ ( J ), one in each of the K intervals. We can compute σ ( J ) directly by noting

1 The spectrum of a Jacobi operator J is given by the closure of the set of E ∈ ❘ for which a nontrivial polynomially bounded solution to J ψ = Eψ exists. When J is periodic, one either has a bounded solution or all solutions grow exponentially on at least one half-line.

that θ = 0 and θ = π give the endpoints of these intervals. We shall thus focus on the two K × K symmetric matrices

$$\i \text{ the two $K\times K$ symmetric matrices} \\ \mathbf J _ { \pm } = \begin{bmatrix} b _ { 1 } & a _ { 1 } & & \pm a _ { K } \\ a _ { 1 } & b _ { 2 } & \ddots & \\ & \ddots & \ddots & \ddots \\ & & \ddots & b _ { K - 1 } & a _ { K - 1 } \\ \pm a _ { K } & & a _ { K - 1 } & b _ { K } \end{bmatrix}. \quad ( 1. 7 ) \\ \text{ecsly, enumerating the $2K\,eigenvalues of $J_{+}$ and $J_{-}$ such that}$$

More precisely, enumerating the 2 K eigenvalues of J + and J -such that

$$E _ { 1 } < E _ { 2 } \leq E _ { 3 } < E _ { 4 } \leq \dots \leq E _ { 2 K - 1 } < E _ { 2 K }$$

(strict inequalities separate eigenvalues from J + and J -), then

$$\sigma ( \mathfrak { J } ) = \bigcup _ { j = 1 } ^ { K } \, [ E _ { 2 j - 1 }, E _ { 2 j } ].$$

This discussion suggests three ways to compute the spectrum (1.9):

- 1. For a fixed E , test if E ∈ σ ( J ) by explicitly calculating tr( M E k ( )) and checking if (1.6) holds;
- 2. Construct the degreeK polynomial tr( M K ( E )) and find the roots of tr( M K ( E )) = ± 2;
- 3. Compute the eigenvalues of the two K × K symmetric matrices J ± .

The first method gives a fast way to test if E ∈ σ ( J ) for a given E , but is an ineffective way to compute the entire spectrum. The second method suffers from the numerical instabilities mentioned earlier. The third approach is most favorable, but O K ( 3 ) complexity and numerical inaccuracies in the computed eigenvalues become a concern for large K .

Jacobi operators with long periods arise as approximations to Schr¨o dinger operators with aperiodic potentials. For one special case (the almost Mathieu potential), Thouless [47] and Lamoureux [29] proposed an O K ( 2 ) algorithm to compute the spectra of J ± for the periodK approximation. Section 2 describes a simpler O K ( 2 ) algorithm that does not exploit special properties of the potential, and so applies to any periodK Jacobi operator. Section 3 addresses aperiodic potentials in some detail, showing how their spectra can be covered by those of periodic approximations. Section 4 shows application of our algorithm to estimate the fractal dimension of the spectrum for aperiodic operators with potential given by primitive substitution rules.

## 2. Algorithm

The conventional algorithm for finding all eigenvalues of a general symmetric matrix A first requires the application of unitary similarity transformations to reduce the matrix to tridiagonal form (in which all entries other than those on the main diagonal and the first super- and sub-diagonals are zero). This transformation is the most costly part of the eigenvalue computation: for a

Figure 1. The nonzero pattern of J ± for K = 32 (left) and an illustration of all the entries that are nonzero at some point in the transformation of J ± to tridiagonal form using the conventional plane rotation approach (right).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-93311d28-56dc-4a34-9c99-4607d302491a.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-2922c2a1-c83a-43cb-a81c-7ddf3c4a123d.jpg)

K × K matrix, the reduction takes O K ( 3 ) operations, while the eigenvalues of the tridiagonal matrix can be found to high precision in O K ( 2 ) further operations [38, § 8.15]. When many entries of A are zero, one might exploit this structure to perform fewer elementary similarity transformations. This is the case when A is banded: a j,k = 0 when | j -k | &gt; b , where b is the bandwidth of A . For fixed b , the tridiagonal reduction takes O K ( 2 ) operations as K →∞ .

The matrix J ± in (1.7) is tridiagonal plus entries in the (1 , K ) and ( K, 1) positions that give it bandwidth b = K . One might still hope to somehow exploit the zero structure to quickly reduce J ± to tridiagonal form. Unfortunately, the transformation that eliminates the (1 , K ) and ( K, 1) entries creates new nonzero entries where once there were zeros, starting a cascade of new nonzeros whose sequential elimination leads again to O K ( 3 ) complexity. Figure 1 shows how the conventional approach to tridiagonalizing a sparse matrix (using Givens plane rotations) eventually creates O K ( 2 ) nonzero entries, requiring O K ( 2 ) storage and O K ( 3 ) floating point operations.

The challenge is magnified by the eigenvalues themselves. As described in Section 3, we approximate aperiodic operators whose spectra are Cantor sets, implying that the eigenvalues of σ ( J ± ) will be tightly clustered. The accuracy of these computed eigenvalues is critical to the applications we envision (e.g., estimating the fractal dimension of the Cantor sets), warranting use of extended (quadruple) precision arithmetic that magnifies the cost of O K ( 3 ) operations. Moreover, these studies often examine a family of operators over many parameter values (e.g., as the { b n } terms are scaled), so expedient algorithms are helpful even when K is tractable for a single matrix.

Thouless [47] and, in later, more complete work, Lamoureux [29] provide an explicit unitary matrix Q ± such that Q J Q ∗ ± ± ± is tridiagonal in the special case of the almost Mathieu potential,

$$a _ { n } = 1, \quad b _ { n } = 2 \lambda \cos ( 2 \pi ( n \alpha + \theta ) ),$$

where α is a rational approximation to an irrational parameter of true interest, and λ and θ are constants. This transformation gives an O K ( 2 ) algorithm for computing σ ( J ± ), but relies on the special structure of the potential.

glyph[negationslash]

Alternatively, we note that one can compute all eigenvalues of J ± in O K ( 2 ) time by simply reordering the rows and columns of J ± to yield a matrix with small bandwidth that is independent of K . Recall that the nonzero pattern of a symmetric matrix A can be represented as an undirected graph having vertices labeled 1 , . . . , K , with distinct vertices j and k joined by an edge if a j,k = a k,j = 0. (We suppress the loops corresponding to a j,j = 0.) The graph for J ± has a single cycle, shown in Figure 2 for K = 8. In the conventional labeling, the corner entries in the (1 , K ) and ( K, 1) positions give an edge between vertices 1 and K . To obtain a matrix with narrow bandwidth, simply relabel the vertices in a breadth-first fashion starting from vertex 1, as illustrated in Figure 2. Now each vertex only connects to vertices whose labels differ by at most two : if we permute the rows and columns of the matrix in accord with the relabeling, the resulting matrix will have bandwidth 2 (i.e., a pentadiagonal matrix). More explicitly, define glyph[negationslash]

$$p ( j ) = \begin{cases} \ 2 j - 1, & j \in \{ 1, \dots, \lceil K / 2 \rceil \} ; \\ 2 ( K - j - 1 ), & j \in \{ \lceil K / 2 \rceil + 1, \dots, K \}. \end{cases}$$

Let P = [ e p (1) , e p (2) , . . . , e p K ( ) ], where e j is the j th column of the K × K identity matrix; then PJ P ± ∗ has bandwidth b = 2. When K = 8,

$$y \, \text{ matrix} ; \, \text{ then } \mathbf P J _ { \pm } \mathbf P ^ { * } \, \text{ has bandwidth } b = 2. \, \text{ When } K = 8, \\ \mathbf J _ { \pm } = \, \left [ \begin{array} { c c c c } b _ { 1 } & a _ { 1 } & & & \pm a _ { 8 } \\ a _ { 1 } & b _ { 2 } & a _ { 2 } & & & \\ & a _ { 2 } & b _ { 3 } & a _ { 3 } & & \\ & a _ { 3 } & b _ { 4 } & a _ { 4 } & & \\ & & a _ { 4 } & b _ { 5 } & a _ { 5 } & \\ & & & a _ { 5 } & b _ { 6 } & a _ { 6 } \\ & & & & a _ { 6 } & b _ { 7 } & a _ { 7 } \\ \pm a _ { 8 } & & & & a _ { 7 } & b _ { 8 } \end{array} \right ] \\ \text{ dered to }$$

is reordered to

$$\text{$s$ reordered to} \\ \text{$PJ_{\pm}P^{*}=\left[\begin{matrix} b _ { 1 } & \pm a _ { 8 } & a _ { 1 } \\ \pm a _ { 8 } & b _ { 8 } & a _ { 7 } \\ a _ { 1 } & b _ { 2 } & a _ { 2 } \\ a _ { 7 } & b _ { 7 } & a _ { 6 } \\ a _ { 2 } & b _ { 3 } & a _ { 3 } \\ a _ { 6 } & b _ { 6 } & a _ { 5 } \\ a _ { 3 } & b _ { 4 } & a _ { 4 } \\ a _ { 5 } & a _ { 4 } & b _ { 5 } \end{matrix} \right ],$$

where unspecified entries are zero. The tridiagonal reduction of banded symmetric matrices has been carefully studied, starting with Rutishauser 2 [42]; see [10] for contemporary algorithmic considerations. Such matrices can be

2 Rutishauser was motivated by pentadiagonal matrices that arise in the addition of continued fractions.

Figure 2. The vertex reordering scheme for K = 8.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-2f11bf03-b4dd-475a-8f6f-95fca8154f06.jpg)

Figure 3. Analogue of Figure 1 for the pentadiagonal matrix from the breadth-first ordering (left) and the nonzeros that arise in reduction of this matrix to tridiagonal form using the bulge-chasing approach (right).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-923675af-81e8-431a-a464-802217dfc08d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-4c8056e2-d87d-4fb4-88cf-af555d753744.jpg)

reduced to tridiagonal form using Givens plane rotations applied in a 'bulgechasing' procedure that increases the bandwidth by one; Figure 3 shows the nonzero entries introduced by this reduction. 3 Removing the ( j + 2 , j ) entry introduces a new entry in the third subdiagonal (a 'bulge') that is 'chased' toward the bottom right with O K ( ) additional Givens rotations, each of which requires O (1) floating point operations. Performing this exercise for j = 1 , . . . , K -2 amounts to O K ( 2 ) work and O K ( ) storage to reduce PJ P ± ∗ to tridiagonal form, an improvement over the usual O K ( 3 ) work and O K ( 2 ) storage. (Our application does not require the eigenvectors of J ± , so we do not store the transformations that tridiagonalize J ± .)

One can compute a banded tridiagonalization using the LAPACK software library's dsbtrd routine, or compute the eigenvalues directly with the

3 An alternative method that finds the eigenvalues of a banded symmetric matrix directly in its band form is described in [38, § 8.16]. This method, based on small Householder reflectors, is particularly effective when a small subset of the spectrum is sought.

Figure 4. Performance of LAPACK's standard ( dsyev ) and banded ( dsbev ) symmetric eigensolvers in double and quadruple precision, applied to the Fibonacci model.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-f65f80ce-c4b2-4787-bcc7-313b1a0d8061.jpg)

banded eigensolver dsbev [2]. (If PJ P ± ∗ stored in sparse format, MATLAB's eig command will identify the band form and find the eigenvalues in O K ( 2 ) time.) We have benchmarked our computations in LAPACK with standard double precision arithmetic and a variant compiled for quadruple precision. 4 When a n = 1 and b n = 0 for all n , the eigenvalues of J ± are known in closed form [23]. With the reordering scheme described above, LAPACK returns the correct eigenvalues, up to the order of machine precision (roughly 10 -16 for double precision and 10 -34 for quadruple precision). Figure 4 compares the timing of the reordered scheme to the traditional dense matrix approach. In both double and quadruple precision, the O K ( 2 ) performance of the reordered approach offers a significant advantage. (These timings were performed on a desktop with a 3.30 GHz Intel Xeon E31245 processor, applied to the Fibonacci model described in the next section. The numbers in the legend indicate the slope of a linear fit of the last five data points, and each algorithm is averaged over the four Fibonacci parameters λ = 1 2 3 4.) , , ,

## 3. Spectral theory for quasiperiodic Schr¨dinger operators o

We focus on Jacobi operators that are discrete Schr¨dinger operators o , that is, those operators for which the off-diagonal terms satisfy a n = 1 for all n , and the potential { b n } varies in a deterministic but non-periodic fashion. 5

4 http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&amp;t=2739 gives compilation instructions.

5 In off-diagonal models, b n = 0 for all n , while the a n coefficients vary aperiodically, see, e.g., [33, 50].

For several prominent examples the spectrum is a Cantor set; in other cases even a gross description of the spectral type has been elusive.

Periodic approximations lead to elegant covers of the spectra of the quasiperiodic operators, 6 and while the arguments that produce these covers are now standard, there is not such a direct venue in the literature where this framework is quickly summarized. Consequently, we recapitulate the essential arguments for the reader unfamiliar with this landscape; those seeking greater detail can consult the more extensive survey [13]. In Section 4 we shall apply our algorithm to numerically compute the spectral covers described here.

Avariety of quasiperiodic potentials have been investigated in the mathematical physics literature; for a survey, see, e.g., [18]. Most prominent is the almost Mathieu operator [25, 47], with potential

$$b _ { n } = 2 \lambda \cos ( 2 \pi ( n \alpha + \theta ) )$$

glyph[negationslash]

for irrational α , nonzero coupling constant λ , and phase θ ∈ ❘ . The spectrum of the almost Mathieu operator is a Cantor set for all irrational α , every θ , and every λ = 0 [3]. Moreover, the Lebesgue measure of the spectrum is precisely 4 1 ∣ ∣ - | λ | ∣ ∣ whenever α is irrational [4, Theorem 1.5]. Of particular interest is the Hausdorff dimension of the spectrum in the critical case λ = 1, about which very little is known; see [30] for some partial results.

Sturmian potentials take the form

$$b _ { n } = \lambda \chi _ { [ 1 - \alpha, 1 ) } ( ( n \alpha \, + \theta ) \bmod 1 ),$$

where χ S is the indicator function on the set S ⊂ ❘ and α , λ , and θ play the same role as in (3.1). For such potentials the spectrum is a Cantor set of zero Lebesgue measure for all irrational α , nonzero λ and phases θ [8]. It is conjectured that for a fixed nonzero coupling λ , the Hausdroff dimension of the spectrum is constant on a set of α 's of full Lebesgue measure; currently, this is only known for λ ≥ 24 [16, Theorem 1.2].

## 3.1. Quasiperiodic operators from substitution rules

Alternatively, aperiodic potentials can be constructed from primitive substitution rules. Unlike the almost Mathieu and Sturmian cases, the spectral type of the operators is not well-understood. Having spectrum of zero Lebesgue measure precludes the presence of absolutely continuous spectrum, a result known well before Damanik and Lenz proved zero-measure spectrum; see the main result of [28]. So far, numerous partial results exclude eigenvalues for particular substitution operators [8, 11, 12, 17, 26], and no results yet establish the existence of eigenvalues for any (two-sided) substitution operator. The construction is slightly more involved than previous examples, but since seeking a better understanding of substitution potentials is a motivation of this work, we describe these objects in detail.

6 Except for the Almost Mathieu operator, the potentials we describe are not quasiperiodic in the classical sense of being almost periodic sequences with finitely generated frequency modules. However, the terminology is completely standard at this point.

↦

↦

↦

↦

↦

↦

↦

Let A ⊂ ❘ be a finite set, called the alphabet . 7 A substitution on A is a rule that replaces elements of A by finite-length words over A . For example, the period doubling substitution, S PD , is defined by the rules a → ab and b → aa over the two-letter alphabet A = { a b , } [9]. The Thue-Morse substitution uses the rules a → ab and b → ba [7, 20]. The Fibonacci potential [27, 36] can be viewed as a Sturmian potential with α = ( √ 5 -1) / 2 and θ = 0, or as a primitive substitution potential with rules a → ab , b → a . (The equivalence of these definitions follows from [8, Lemma 1b].) One can also study substitutive sequences on larger alphabets; for example, the Rudin-Shapiro substitution is defined on the four-symbol alphabet A = { a b c d , , , } by the rules a → ab , b → ac , c → db , and d → dc [41, 43].

↦

↦

Each of the aforementioned examples enjoys a property known as primitivity, which can be described informally as the existence of an iterate which maps each letter to a word containing the full alphabet. More precisely, we say that a substitution S is primitive if there exists some k ∈ ❩ + so that for every a b , ∈ A , b is a subword of S k ( a ).

↦

We now describe how a primitive substitution rule leads to a quasiperiodic Schr¨dinger operator, focusing on period doubling as a concrete example. o To obtain an aperiodic sequence, start with the symbol a and form the sequence w k = S k PD ( a ) by iteratively applying the period doubling substitution rules, a → ab and b → aa . The result is a sequence of finite words: w 0 = a , w 1 = ab , w 2 = abaa , w 3 = abaaabab , and so on. Notice that w k is always a prefix of w k +1 . In particular, there is a well-defined limiting sequence,

↦

$$x _ { P D } = \lim _ { k \to \infty } S _ { P D } ^ { k } ( a ) = a b a a b a b a b a a a a \dots, \quad \quad ( 3. 3 )$$

with the property that x PD is fixed by the period doubling substitution; sequences with this property are called substitution words for the substitution.

From a substitution word x PD one can construct a quasiperiodic potential for a discrete Schr¨dinger operator. Given a coupling constant o λ , take

$$a = \lambda ; \quad b = 0 ; \quad b _ { n } = n t h \, \text{symbol of } x _ { P D }. \quad \quad ( 3. 4 )$$

One subtlety remains: the substitution word x PD is a one-sided sequence, while our potentials are two-sided, i.e., we should specify b n for all n ∈ ❩ . To generate two-sided potentials, one considers an accumulation point of leftshifts of x PD ; equivalently, one considers a two-sided sequence with the same local factor structure as x PD . The details of this construction follow.

Suppose S is a primitive substitution on the finite alphabet A ⊂ ❘ and x ∈ A ◆ is some substitution word thereof, i.e., S x ( ) = x . The subshift generated by S is the set of all two-sided sequences with the same local factor structure as x . More precisely,

Ω S = { ω ∈ A ❩ : ω n · · · ω m is a subword of x for every n ≤ m . } (3.5) Using primitivity of S , one can check that this definition of Ω S does not depend upon the choice of substitution word. One can also generate the set

7 We do not necessarily have to restrict our attention to alphabets consisting of real numbers, but this makes the definition of our subshift potentials in (3.7) somewhat simpler.

↦

Ω S via the following dynamical procedure. First, endow the sequence space A ❩ with some metric that induces the product topology thereupon, e.g.,

❩ where δ a,b denotes the usual Kronecker delta symbol. The sequences ω and ω ′ are close with respect to the metric in (3.6) if and only if they agree on a large window centered at the origin, so this metric does indeed induce the topology of pointwise convergence on Ω S . One can check that Ω S defined by (3.5) coincides precisely with the set of limits of convergent subsequences of the sequence ( T n x ) ∞ n =1 , where T denotes the usual left shift ( Tω ) n = ω n +1 . (There is a minor technicality here: since x is a one-sided sequence ( T n x ) k is only well-defined for n &gt; | k | when k is negative.)

$$d ( \omega, \omega ^ { \prime } ) = \sum _ { n \in \mathbb { Z } } \frac { 1 - \delta _ { \omega _ { n }, \omega _ { n } ^ { \prime } } } { 2 ^ { | n | + 1 } },$$

Primitivity of S implies that the topological dynamical system (Ω S , T ) is strictly ergodic [39, Proposition 5.2 and Theorem 5.6]. For each ω ∈ Ω , S one obtains a Schr¨dinger operator o H ω on glyph[lscript] 2 ( ❩ ) defined via

$$( H _ { \omega } \psi ) _ { n } = \psi _ { n - 1 } + \omega _ { n } \psi _ { n } + \psi _ { n + 1 }$$

for each n ∈ ❩ . (The right hand side of (3.7) makes sense, since we chose A ⊆ ❘ .) Minimality of (Ω S , T ) implies ω -invariance of the spectrum.

Proposition 3.1. Given Ω S and H ω as above, there is a uniform compact set Σ with the property that σ H ( ω ) = Σ for every ω ∈ Ω S .

Proof. Given ω, ω ′ ∈ Ω , there is a sequence S { n k } ⊂ ❩ with the property that T n k ω → ω ′ as k →∞ . (This is a consequence of minimality of (Ω S , T ), which follows from [39, Proposition 4.7], for example.) In particular,

$$H _ { \omega ^ { \prime } } = s _ { \substack { k \rightarrow \infty \\ k \rightarrow \infty } } H _ { T ^ { n } _ { k } \omega }.$$

By a standard strong approximation argument (e.g., [40, Theorem VIII.24]),

$$\sigma ( H _ { \omega ^ { \prime } } ) \subseteq \bigcap _ { \ell = 1 } ^ { \infty } \overline { \bigcup _ { k = \ell } ^ { \infty } \sigma ( H _ { T ^ { n _ { k } } \omega } ) } = \sigma ( H _ { \omega } ).$$

By symmetry, we can run the previous argument with the roles of ω and ω ′ reversed, so σ H ( ω ) = σ H ( ω ′ ). glyph[square]

This is the upshot of the previous proposition: if we want to study the spectrum of a two-sided substitutive Schr¨dinger operator o as a set , then we can work with any member of the associated subshift.

One can avoid the construction needed to generate two-sided potentials by considering instead half-line Jacobi operators, J + : glyph[lscript] 2 ( ◆ ) → glyph[lscript] 2 ( ◆ ). These operators are defined entrywise by (1.1) for n ∈ ◆ with the normalization ψ 0 = 0 to make ( J + ψ ) 1 well defined. If J + denotes the Jacobi operator defined by a n ≡ 1 and (3.4), and J denotes a two-sided Jacobi operator generated in one of the two fashions just described, then by [44, Theorem 7.2.1],

$$\sigma _ { \text{ess} } ( \mathfrak { I } _ { + } ) = \sigma ( \mathfrak { I } ),$$

where σ ess ( ) · denotes the essential spectrum. In particular, σ ( J + ) and σ ( J ) have the same Hausdorff dimension, since (3.9) implies that they coincide up to a countable set of isolated points.

## 3.2. Periodic approximations

All the classes of quasiperiodic models we have described have natural periodic approximations. For the almost Mathieu and Sturmian cases: replace the irrational α with a rational approximant; for the substitution rules: pick a starting symbol, generate a string from finitely many applications of the substitution rules, and repeat that string periodically. We seek high fidelity numerical approximations to these periodic spectra, as a vehicle for understanding properties of quasiperiodic models, such as the fractal dimension of the spectrum.

In the case of substitution rules, we shall describe how periodic approximations lead to an upper bound on the spectrum of the quasiperiodic operator. Again we focus on the period doubling substitution. Fix S PD , x PD as before, choose some ω ∈ Ω PD , and put w a k = S k ( a ) and w b k = S k ( b ). We generate periodic Schr¨dinger operators o H a k and H b k by repeating the strings w a k and w b k periodically, giving H a k and H b k with coefficients having period 2 k . Moreover, we choose H a k and H b k in such a way that

$$H _ { \omega } = \underset { k \rightarrow \infty } { s \text{-} \lim } H _ { k } ^ { a } = \underset { k \rightarrow \infty } { s \text{-} \lim } H _ { k } ^ { \text{b} }.$$

Define Σ a k ≡ σ H ( a k ) and Σ b k ≡ σ H ( b k ) for all k ≥ 0. Application of strong approximation, as before, implies

$$\Sigma \subseteq \bigcap _ { n = 1 } ^ { \infty } \overline { \bigcup _ { k = n } ^ { \infty } \Sigma _ { k } ^ { a } }, \quad \Sigma \subseteq \bigcap _ { n = 1 } ^ { \infty } \overline { \bigcup _ { k = n } ^ { \infty } \Sigma _ { k } ^ { b } }.$$

These unions are too large to be computationally tractable, but the hierarchical structure of the periodic approximations saves the day. In these cases the monodromy matrices (1.3) take the forms

$$M _ { k } ^ { a } ( E ) \ & = \ \begin{bmatrix} E - w _ { k } ^ { a } ( 2 ^ { k } ) & - 1 \\ 1 & 0 \end{bmatrix} \cdots \begin{bmatrix} E - w _ { k } ^ { a } ( 1 ) & - 1 \\ 1 & 0 \end{bmatrix} \\ M _ { k } ^ { b } ( E ) \ & = \ \begin{bmatrix} E - w _ { k } ^ { b } ( 2 ^ { k } ) & - 1 \\ 1 & 0 \end{bmatrix} \cdots \begin{bmatrix} E - w _ { k } ^ { b } ( 1 ) & - 1 \\ 1 & 0 \end{bmatrix},$$

$$\dots$$

where now the subscript k denotes the k th iteration of the substitution rule, and hence a full period of length K = 2 . Recalling the discriminant condik tion (1.6), define

$$x _ { k } ( E ) \equiv \text{tr} ( M _ { k } ^ { a } ( E ) ), \ \ y _ { k } ( E ) \equiv \text{tr} ( M _ { k } ^ { b } ( E ) ).$$

As in (1.6), these functions encode the spectra of the periodic approximants, in that Σ a k = { E : | x k ( E ) | ≤ } 2 and Σ b k = { E : | y k ( E ) | ≤ } 2 . The rules of the

period doubling substitution imply

$$M _ { k + 1 } ^ { a } = M _ { k } ^ { \text{b} } M _ { k } ^ { a }$$

$$M _ { k + 1 } ^ { \flat } = M _ { k } ^ { a } M _ { k } ^ { a }.$$

Applying the Cayley-Hamilton theorem ( M 2 k -tr( M M k ) k + I = 0) to (3.11) and (3.12) yields

$$x _ { k + 1 } = x _ { k } y _ { k } - 2$$

$$y _ { k + 1 } = x _ { k } ^ { 2 } - 2,$$

given in [9, eq. (1.9)]. Now (3.13) and (3.14) imply Σ a k +1 ∪ Σ b k +1 ⊆ Σ a k ∪ Σ b k for all k ≥ 0, thus reducing (3.10) to the more tractable

$$[ \text{Period Doubling} ] \quad \ \ \Sigma \subseteq \left \{ \, \, \right \} \Sigma _ { k } ^ { a } \cup \Sigma _ { k } ^ { b }.$$

$$\Sigma \subseteq \bigcap _ { k = 1 } ^ { \infty } \Sigma _ { k } ^ { a } \cup \Sigma _ { k } ^ { b }.$$

The spectrum for the Thue-Morse substitution can be expressed in the same way, except that substitution rule replaces (3.13)-(3.14) with

$$x _ { k + 1 } = x _ { k - 1 } ^ { 2 } ( x _ { k } - 2 ) + 2 \quad \quad \quad$$

$$y _ { k + 1 } & = x _ { k + 1 }.$$

Reasoning as in the period doubling case, (3.16)-(3.17) implies

[Thue-Morse]

$$\Sigma \subseteq \bigcap _ { k = 1 } ^ { \infty } \Sigma _ { k } ^ { a } \cup \Sigma$$

Notice that (3.15) and (3.18) are not identical. In general, the covers of the spectrum obtained via this approach depend quite strongly on the chosen substitution. Note, however, that the same formula (3.18) holds for the Fibonacci substitution. The recursive relationships in (3.13)-(3.14) and (3.16)(3.17) are known as the trace maps for the period doubling and Thue-Morse potentials. Similar trace maps can be constructed for arbitrary substitutions; see [5, Theorem 1] and [6, Theorem 1] for their description.

Figure 5 depicts the covers Σ a k ∪ Σ b k for the period doubling potential and Σ a k ∪ Σ a k +1 for Thue-Morse as λ increases from zero with k = 7, giving an impression of how rapidly the covering intervals shrink as λ increases. Notice the quite different nature of the spectra obtained from these two substitution rules. Similar figures for the Fibonacci spectrum are shown in [14]. (The algorithm described in the last section enables such computations for much larger values of k , but the resulting figures become increasingly difficult to render due to the number and narrow width of the covering intervals.)

## 4. Numerical calculation of periodic covers

In this section, we apply the algorithm described in Section 2 to compute the spectral covers for quasiperiodic operators from substitution rules described in Section 3. We begin with some calculations that emphasize the need for

Figure 5. Spectral covers for the period doubling (Σ a k ∪ Σ ) b k and Thue-Morse (Σ a k ∪ Σ a k +1 ) potentials as a function of the coupling constant λ , for k = 7.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-9aa04470-07a3-42b4-9f29-48114ea80171.jpg)

higher precision arithmetic to resolve the spectrum, then estimate various spectral quantities for these quasiperiodic operators.

## 4.1. Necessity for extended precision

High fidelity spectral approximations for quasiperiodic Schr¨dinger operators o are tricky to compute due not only to the complexity of the eigenvalue problem, but also considerations of numerical precision. The LAPACK symmetric eigenvalue algorithms are expected to compute eigenvalues of the K × K matrix J ± accurate to within p K ( ) ‖ J ± ‖ ε mach [2, p. 104], where p K ( ) is a 'modestly growing function of K ' and ε mach denotes the machine epsilon value for the floating point arithmetic system (on the order of 10 -16 for double precision and 10 -34 for quadruple precision [1]). Figure 6 shows the width of the smallest interval in the approximation Σ a k for the period doubling and Thue-Morse potentials, as computed in double and quadruple precision arithmetic for various values of λ . As λ increases, this smallest interval shrinks ever quicker. Comparing double and quadruple precision values, one sees that for even moderate values of K = 2 , double precision is unable to accurately k

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-3a3ab8d6-829e-4e16-9109-3991b1e675c4.jpg)

k

Figure 6. Minimum interval length in Σ a k for a range of λ values (powers of two) for the period doubling and ThueMorse potentials. The solid line shows double precision, the dashed line quadruple precision. Black data points are believed to be correct to plotting accuracy; gray ones are not.

resolve the spectrum. 8 (Similar numerical errors will be observed for the original matrix J ± . Since it only permutes matrix entries, our algorithm does not introduce any new instabilities.)

8 The calculations of fractal dimension and gaps that follow were performed in quadruple precision arithmetic, and we generally restricted the values of λ and k to obtain numerically reliable results. In the event the numerical results violate the ordering of eigenvalues of J + and J -in (1.8), any offending interval is replaced by one of width 20 ε mach centered at the midpoint of the computed eigenvalues.

## 4.2. Approximating fractal dimensions

The quasiperiodic spectrum Σ is known to be a Cantor set for the Fibonacci, period doubling, and Thue-Morse potentials for all λ &gt; 0, suggesting calculation of the fractal dimensions of these sets as functions of λ . We begin with two standard definitions; see, e.g., [22, 37].

Definition 4.1. Given A ⊂ ❘ and some α ∈ [0 , 1], let

$$h ^ { \alpha } ( A ) \equiv \lim _ { \Delta \to 0 } \inf _ { \Delta \text{-covers} } \sum$$

where a ∆-cover of A is defined to be collection { B m m } ≥ 1 of intervals such that A ⊂ ⋃ m ≥ 1 B m and | B m | &lt; ∆ for each m . The Hausdorff dimension of A is then

$$\dim _ { H } ( A ) \equiv \inf \{ \alpha \, \colon h ^ { \alpha } ( A ) < \infty \}$$

Definition 4.2. Given A ⊂ ❘ , let N A ( ε ) denote the number of intervals of the form [ jε, ( j +1) ), ε j ∈ ❩ , that have a nontrivial intersection with A . The upper and lower box counting dimensions of A are

$$\dim _ { B } ^ { \pm } ( A ) \equiv \lim _ { \varepsilon \to 0 } \sup _ { \inf } \ \frac { \log N _ { A } ( \varepsilon ) } { \log 1 / \varepsilon } ;$$

when these agree, the result is the box-counting dimension of A , dim B ( A ).

The analysis in the last section provides natural covering sets for Σ, which we denote as C k (i.e., for period doubling, C k ≡ Σ a k ∪ Σ ; for Fibonacci b k and Thue-Morse, C k ≡ Σ a k ∪ Σ a k +1 ). We estimate dim H (Σ) using a heuristic proposed by Halsey et al. [24]. Enumerate the n k intervals in C k as { B k,m } :

$$C _ { k } = \bigcup _ { m = 1 } ^ { n _ { k } } B _ { k, m }.$$

## Approximating the Hausdorff dimension

- 1: Input : two consecutive covers C k,λ , C k +1 ,λ of the Cantor spectrum.
- 2: Construct the function f ( α ) ≡ ∑ n k m =1 | B k,m | α -∑ n k +1 m =1 | B k +1 ,m | α .
- 3: Compute root α k of f ( α ) in the interval [0 , 1].
- 4: Output : α k as the approximation to the dimension.

For the Fibonacci case, Damanik et al. [15] proved upper and lower bounds on dim H (Σ) for λ ≥ 8 in terms of the functions S u ( λ ) = 2 λ +22 and S l ( λ ) = 1 2 ( λ -4 + √ ( λ -4) 2 -12 . To benchmark our method, we compare ) the approximate dimension with the upper and lower bounds, as shown in Figure 7, indeed obtaining satisfactory results. Figure 8 shows estimates to dim H (Σ) drawn from spectral covers for k = 6 , . . . , 10, with good convergence in k in the small λ regime. The narrow covering intervals for large λ and k pose a significant numerical challenge, even in quadruple precision arithmetic.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-60212ba0-dcb7-4987-94eb-b6da80576828.jpg)

λ

Figure 7. Estimates of dim H (Σ) for the Fibonacci operator (solid line) computed from Σ a k ∪ Σ a k +1 ( k = 15 16), , obeying the upper and lower bounds from [15].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-38447e27-5aac-42c0-a077-4bcb85a49df4.jpg)

λ

Figure 8. Estimates of dim H (Σ) for the period-doubling operator using consecutive covers Σ a k ∪ Σ b k and Σ a k +1 ∪ Σ b k +1 , for k = 6 , . . . , 10. (The k = 9 and k = 10 plots do not show the full range of λ values, due to the numerical challenge of working with large k and λ .)

The bounds from [15] imply dim H (Σ) → 0 like log(1 + √ 2) / log( λ ) as λ →∞ for the Fibonacci case. In contrast, Liu and Qu [32] recently showed that for the Thue-Morse operator, dim H (Σ) is bounded away from zero as λ → ∞ . Interestingly, the approximation scheme described above does not yield consistent results for Thue-Morse, perhaps a reflection of the exotic behavior identified in [32]. As an alternative, in Figure 9 we show estimates of the box-counting dimension for Thue-Morse. Since for any finite k the covers comprise the finite union of real intervals, log( N C k ( ε )) / log(1 /ε ) → 1 as ε → 0. However, for finite ε the shape of the curves in Figure 9 can suggest rough values for dim B (Σ); cf. [45]. Figure 10 repeats the experiment for period doubling; note that the plot for λ = 4 shows good agreement with the estimate for dim H (Σ) seen in Figure 8. (The Hausdorff and box counting

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-085220c3-201c-4e53-9ebb-ff302026f35d.jpg)

λ

ε

Figure 9. Estimates of dim B (Σ) for the Thue-Morse operator for λ = 4 8, using covers , C k with k = 5 , . . . , 9.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-bdb94cfe-f64f-428f-a36a-9f8f85b138b4.jpg)

Figure 10. Estimates of dim B (Σ) for the period-doubling operator for λ = 4 8, using covers , C k with k = 5 , . . . , 9.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-c15fdec0-8204-4cba-a2dd-09c8de97be84.jpg)

dimensions need not be equal; considerable effort was devoted over the years to showing dim H (Σ) = dim B (Σ) for all λ &gt; 0 for the Fibonacci case, with the complete result obtained only recently [19].)

The spectral covers we have described can behave in rather subtle ways. To illustrate this fact, Figure 11 shows the largest gap in the Thue-Morse cover Σ a k ∪ Σ a k +1 as a function of the parameter λ . Bellisard [7] showed that, for the aperiodic model, this gap should behave like λ log 4 / log 3 as λ → 0. The covers satisfy this characterization for intermediate values of λ , but appear to behave instead like λ 2 as λ → 0. The larger the value of k (hence the longer the periodic approximations), the larger the range of λ values that are consistent with Bellisard's spectral asymptotics. This scenario provides another justification for the study of largek approximations: for this potential, one must consider large k to approximate the spectrum adequately for small λ .

Our algorithm also expedites study of the square Hamiltonians J sq acting on glyph[lscript] 2 ( ❩ ) × glyph[lscript] 2 ( ❩ ) via

$$( \mathfrak { I } _ { s q } \psi ) _ { m, n } = \psi _ { m - 1, n } + \psi _ { m, n - 1 } + ( b _ { m } + b _ { n } ) \psi _ { m, n } + \psi _ { m + 1, n } + \psi _ { m, n + 1 }.$$

= 8

Figure 11. Computed values of the largest gap in the cover Σ a k ∪ Σ a k +1 of the Thue-Morse spectrum, along with Bellisard's asymptotic description λ log 4 / log 3 . As k increases, the covers obey the asymptotics for smaller values of λ .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-8304b501-8be0-4290-b730-cc673b3469c1.jpg)

The spectrum of this operator is Σ + Σ, where Σ is the spectrum of the standard 1-dimensional operator with potential { b n } . When Σ is a Cantor set, Σ + Σ could be an interval, a union of intervals, a Cantorval , or a Cantor set; see, e.g., [34, 35, 37]. For further details and examples for the Fibonacci case, see [14] and the computations, based on σ ( J ± ), of Even-Dar Mandel and Lifshitz [21]. For the period doubling and Thue-Morse potentials, Figure 12 illustrates how the covers of the spectra for these square Hamiltonians develop in λ , revealing values of λ for which Σ + Σ cannot comprise a single interval.

## 5. Conclusion

We have presented a simple O K ( 2 ) algorithm to numerically compute the spectrum of a periodK Jacobi operator that can be implemented in a few lines of code, then we used this algorithm to estimate spectral quantities associated with quasiperiodic Schr¨dinger operators derived from substituo tion rules. The algorithm facilitates study of long-period approximations and extensive parameter studies across a family of operators. Yet even with the efficient algorithm, the accuracy of the numerically computed eigenvalues must be carefully monitored, or even enhanced with extended precision computations. This tool can expedite numerical experiments to help formulate conjectures about spectral properties for a range of quasiperiodic operators.

## Acknowledgments

We thank a referee for helpful comments that led us to study the largest Thue-Morse gap (Figure 11), and David Damanik, Paul Munger, Beresford Parlett, and Dan Sorensen for fruitful discussions. We especially thank May Mei for suggesting computations with substitution operators and sharing numerical results, and Adrian Forster for his Thue-Morse computations.

Figure 12. Spectral covers for the square Hamiltonian for period doubling and Thue-Morse substitutions, with k = 4.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0370v2-figure-0438b510-0927-4360-98f0-14f878d9736c.jpg)

## References

- [1] IEEE Standard for Floating-Point Arithmetic (IEEE Standard 754-2008). Institute of Electrical and Electronics Engineers, Inc., 2008.
- [2] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK User's Guide . SIAM, Philadelphia, third edition, 1999.
- [3] A. Avila and S. Jitomirskaya. The ten martini problem. Ann. of Math. , 170:303-342, 2009.
- [4] A. Avila and R. Krikorian. Reducibility or nonuniform hyperbolicity for quasiperiodic Schr¨dinger cocycles. o Ann. of Math. , 164:911-940, 2006.
- [5] Y. Avishai and D. Berend. Trace maps for arbitrary substitution sequences. J. Phys. A , 26:2437-2443, 1993.
- [6] Y. Avishai, D. Berend, and D. Glaubman. Minimum-dimension trace maps for substitution sequences. Phys. Rev. Lett. , 72:1842-1845, 1994.
- [7] J. Bellissard. Spectral properties of Schr¨dinger's operator with a Thue-Morse o potential. In Number Theory and Physics , volume 47 of Springer Proceedings in Physics , pages 140-150. Springer-Verlag, Berlin, 1990.

- [8] J. Bellissard, B. Iochum, E. Scoppola, and D. Testard. Spectral properties of one-dimensional quasicrystals. Comm. Math. Phys. , 125:527-543, 1989.
- [9] Jean Bellissard, Anton Bovier, and Jean-Michel Ghez. Spectral properties of a tight binding Hamiltonian with period doubling potential. Comm. Math. Phys. , 135:379-399, 1991.
- [10] Christian H. Bischof, Bruno Lang, and Xiaobai Sun. A framework for symmetric band reduction. ACM Trans. Math. Software , 26:581-601, 2000.
- [11] D. Damanik. Singular continuous spectrum for a class of substitution Hamiltonians II. Lett. Math. Phys. , 54:25-31, 2000.
- [12] D. Damanik. Uniform singular continuous spectrum for the period doubling Hamiltonian. Ann. Henri Poincar´ e , 2:101-118, 2001.
- [13] D. Damanik. Strictly ergodic subshifts and associated operators. In Spectral Theory and Mathematical Physics: a Festschrift in Honor of Barry Simon's 60th Birthday , volume 76 of Proc. Sympos. Pure Math. , pages 539-563. American Mathematical Society, Providence, RI, 2007.
- [14] D. Damanik, M. Embree, and A. Gorodetski. Spectral properties of Schr¨dinger o operators arising in the study of quasicrystals, 2012. arXiv:1210.5753 [math.SP].
- [15] D. Damanik, M. Embree, A. Gorodetski, and S. Tcheremchantsev. The fractal dimension of the spectrum of the Fibonacci Hamiltonian. Comm. Math. Phys. , 280:499-516, 2008.
- [16] D. Damanik and A. Gorodetski. Almost sure frequency independence of the dimension of the spectrum of Sturmian Hamiltonians, 2014. arXiv:1406.4810 [math.SP].
- [17] D. Damanik and D. Lenz. Uniform spectral properties of one-dimensional quasicrystals, I. Absence of eigenvalues. Comm. Math. Phys. , 207:687-696, 1999.
- [18] David Damanik and Jake Fillman. Spectral Theory of Discrete OneDimensional Ergodic Schr¨dinger Operators o . In preparation.
- [19] David Damanik, Anton Gorodetski, and William Yessen. The Fibonacci Hamiltonian, 2014. arXiv:1403.7823 [math.SP].
- [20] F. Delyon and J. Peyri`re. Recurrence of the eigenstates of a Schr¨dinger ope o erator with automatic potential. J. Stat. Phys. , 64:363-368, 1991.
- [21] S. Even-Dar Mandel and R. Lifshitz. Electronic energy spectra of square and cubic Fibonacci quasicrystals. Phil. Mag. , 88:2261-2273, 2008.
- [22] Kenneth Falconer. Fractal Geometry: Mathematical Foundations and Applications . Wiley, Chichester, 2nd edition, 2003.
- [23] C. W. Gear. A simple set of test matrices for eigenvalue programs. Math. Comp. , 23:119-125, 1969.
- [24] Thomas C. Halsey, Mogens H. Jensen, Leo P. Kadanoff, Itamar Procaccia, and Boris I. Shraiman. Fractal measures and their singularities: the characterization of strange sets. Phys. Rev. A , 33:1141-1151, 1986.
- [25] P. G. Harper. Single band motion of conduction electrons in a uniform magnetic field. Proc. Phys. Soc. A , 68:874-878, 1955.
- [26] A. Hof, O. Knill, and B. Simon. Singular continuous spectrum for palindromic Schr¨dinger operators. o Comm. Math. Phys. , 174:149-159, 1995.

- [27] Mahito Kohmoto, Leo P. Kadanoff, and Chao Tang. Localization problem in one dimension: mapping and escape. Phys. Rev. Lett. , 50:1870-1872, 1983.
- [28] S. Kotani. Jacobi matrices with random potentials taking finitely many values. Rev. Math. Phys. , 1:129-133, 1989.
- [29] Michael P. Lamoureux. Reflections on the almost Mathieu operator. Integral Equations Operator Theory , 28:45-59, 1997.
- [30] Y. Last. Zero measure spectrum for the almost Mathieu operator. Comm. Math. Phys. , 164:421-432, 1994.
- [31] Yoram Last. Quantum dynamics and decompositions of singular continuous spectra. J. Funct. Anal. , 142:406-445, 1996.
- [32] Qinghui Liu and Yanhui Qu. Iteration of polynomial pair under Thue-Morse dynamic, 2014. arXiv:1403.2257 [math.DS].
- [33] Laurent Marin. On- and off-diagonal Sturmian operators: Dynamic and spectral dimension. Rev. Math. Phys. , 24(05):1250011, 2012.
- [34] Pedro Mendes and Fernando Oliveira. On the topological structure of the arithmetic sum of two Cantor sets. Nonlinearity , 7:329-343, 1994.
- [35] Carlos G. Moreira, Eduardo Mu˜oz Morales, and Juan Rivera-Letelier. On the n topology of arithmetic sums of regular Cantor sets. Nonlinearity , 13:2077-2087, 2000.
- [36] Stellan Ostlund, Rahul Pandit, David Rand, Hans Joachim Schellnhuber, and Eric D. Siggia. One-dimensional Schr¨dinger equation with an almost periodic o potential. Phys. Rev. Lett. , 50:1873-1876, 1983.
- [37] Jacob Palis and Floris Takens. Hyperbolicity and Sensitive Chaotic Dynamics at Homoclinic Bifurcations . Cambridge University Press, Cambridge, 1993.
- [38] Beresford N. Parlett. The Symmetric Eigenvalue Problem . SIAM, Philadelphia, SIAM Classics edition, 1998.
- [39] M. Queff´ elec. Substitution Dynamical Systems - Spectral Analysis . Springer, Berlin, 1987.
- [40] Michael Reed and Barry Simon. Methods of Modern Mathematical Physics I: Functional Analysis . Academic Press, San Diego, revised and enlarged edition, 1980.
- [41] W. Rudin. Some theorems on Fourier coefficients. Proc. Amer. Math. Soc. , 10:855-859, 1959.
- [42] H. Rutishauser. On Jacobi rotation patterns. In Experimental Arithmetic, High Speed Computing and Mathematics , volume 15 of Proceedings of Symposia in Applied Mathematics , pages 219-239. American Mathematical Society, Providence, RI, 1963.
- [43] H. S. Shapiro. Extremal problems for polynomials and power series. Master's thesis, Massachusetts Institute of Technology, 1951.
- [44] Barry Simon. Szeg˝'s o Theorem and Its Descendants . Princeton University Press, Princeton, NJ, 2011.
- [45] Tam´ as T´l, Agnes F¨l¨p, and Tam´s Vicsek. Determination of fractal dimene ´ u o a sions for geometrical multifractals. Physica A , 159:155-166, 1989.
- [46] Gerald Teschl. Jacobi Operators and Completely Integrable Nonlinear Lattices . American Mathematical Society, Providence, RI, 2000.

- [47] D. J. Thouless. Bandwidths for a quasiperiodic tight-binding model. Phys. Rev. B , 28:4272-4276, 1983.
- [48] Morikazu Toda. Theory of Nonlinear Lattices . Springer-Verlag, Berlin, 2nd edition, 1989.
- [49] Pierre van Moerbeke. The spectrum of Jacobi matrices. Inv. Math. , 37:45-81, 1976.
- [50] W.N. Yessen. Spectral analysis of tridiagonal Fibonacci Hamiltonians. J. Spectral Theory , 3:101-128, 2013.

## Charles Puelz

Department of Computational and Applied Mathematics

Rice University

6100 Main Street - MS 134

Houston, Texas

77005-1892

USA

e-mail:

cpuelz@rice.edu

## Mark Embree

Department of Mathematics

Virginia Tech

225 Stanger Street - 0123

Blacksburg, Virginia

24061

USA

e-mail:

embree@vt.edu

Jake Fillman Department of Mathematics Rice University 6100 Main Street - MS 136 Houston, Texas 77005-1892 USA e-mail: jdf3@rice.edu
Charles Puelz Mark Embree Jake Fillman
Spectral Approximation for Quasiperiodic Jacobi Operators
Quasiperiodic Jacobi operators arise as mathematical models of quasicrystals
and in more general studies of structures exhibiting aperiodic order. The
spectra of these self-adjoint operators can be quite exotic, such as Cantor
sets, and their fine properties yield insight into associated dynamical
systems. Quasiperiodic operators can be approximated by periodic ones, the
spectra of which can be computed via two finite dimensional eigenvalue
problems. Since long periods are necessary to get detailed approximations, both
computational efficiency and numerical accuracy become a concern. We describe a
simple method for numerically computing the spectrum of a period-$K$ Jacobi
operator in $O(K^2)$ operations, and use it to investigate the spectra of
Schr\"odinger operators with Fibonacci, period doubling, and Thue-Morse
potentials.