Universal algorithms for computing spectra of periodic operators

Schrödinger operators with periodic (possibly complex-valued) potentials and discrete periodic operators (possibly with complex-valued entries) are considered, and in both cases the computational spectral problem is investigated: namely, under what conditions can a ‘one-size-fits-all’ algorithm for computing their spectra be devised? It is shown that for periodic banded matrices this can be done, as well as for Schrödinger operators with periodic potentials that are sufficiently smooth. In both cases implementable algorithms are provided, along with examples. For certain Schrödinger operators whose potentials may diverge at a single point (but are otherwise well-behaved) it is shown that there does not exist such an algorithm, though it is shown that the computation is possible if one allows for two successive limits.


Introduction and Main Results
We study the computational spectral problem for periodic discrete operators, acting in 2 (Z), as well as Schrödinger operators with periodic potentials acting in L 2 (R d ).We show that it is possible to compute their respective spectra as limits of finite-dimensional approximations.However, in the Schrödinger case this becomes impossible if the potential is allowed to be discontinuous at a single point (but otherwise it is smooth).More precisely, we prove: Theorem 1.1 (Periodic banded matrices).Let B( 2 (Z)) denote the set of bounded operators on 2 (Z) and let {e i } i∈Z be a basis for 2 (Z).Let Ω per ⊂ B( 2 (Z)) be the class of banded, periodic operators with respect to the basis {e i } i∈Z and Ω per N,b the subset of matrices with period N and bandwidth b.Then (i) there exists an algorithm that can compute the spectrum σ(A) of any A ∈ Ω per as the limit of a sequence computable approximations; (ii) there exists an algorithm that can compute the spectrum σ(A) of any A ∈ Ω per N,b with guaranteed error bounds.
Theorem 1.2 (Schrödinger: good case).For d ∈ N define the class of potentials Ω Sch := {V : R d → C | V is 1-periodic and V | (0,1) d ∈ W 1,p ((0, 1) d ) for some p > d}, and given M > 0 and p > d define the class Ω Sch p,M := {V ∈ Ω Sch | V W 1,p ((0,1) d ) ≤ M }.Then (i) there exists an algorithm that can compute the spectrum σ(H) of any Schrödinger operator H = −∆ + V with V ∈ Ω Sch as the limit of a sequence computable approximations; (ii) moreover, for V ∈ Ω Sch p,M this algorithm yields spectral inclusion with guaranteed error bounds.
Theorem 1.3 (Schrödinger: bad case).Let x 0 ∈ [0, 1] and let \ {x 0 } .Then there does not exist an algorithm that can compute the spectrum σ(H) of any Schrödinger operator H = −∆+V with V ∈ Ω Sch x0 as the limit of a sequence computable approximations.However, there does exist an algorithm that can compute σ(H) by taking two successive limits.
Below, in Section 2, we give precise definitions of what an 'algorithm' is,what information is available to it, how it computes, and hence what we mean when we say 'computable'.Informally, Theorems 1.1 and 1.2 imply that these computations can be performed numerically.In fact, we provide actual algorithms which can access the matrix entries (in Theorem 1.1) and pointwise evaluations of the potential (in Theorem 1.2).This should be contrasted with Theorem 1.3 which implies that it is impossible to devise an algorithm that can compute the spectrum of any Schrödinger operator with a potential belonging to Ω Sch x0 .We prove this by contradiction: assuming the existence of such an algorithm, we explicitly construct a potential V ∈ Ω Sch x0 for which this algorithm would fail in its attempt to compute the spectrum.
These statements are nontrivial.The existence of a 'one-size-fits-all' algorithm as in Theorems 1.1 and 1.2 is not obvious: while there are techniques for computing the spectrum of a given operator in each of the above cases -these are indeed well-studied problems -here we prove the existence of a single algorithm (which can be coded) that can handle any input from the given class, without any additional a priori information.Moreover, in the cases of Ω per N,b and Ω Sch p,M there are even guaranteed error bounds.Conversely, the non-existence result of Theorem 1.3 is also not obvious: we prove that regardless of what operations are allowed, as long as the algorithm can only read a finite amount of information at each iteration, there will necessarily be a potential for which the computation will fail.In the spirit of the Solvability Complexity Theory (see a brief discussion below in Subsection 1.1, with more details in Section 2) we show that if one allows for two successive limits (which cannot be collapsed to a single limit), the computation is possible.Note that although the class of potentials Ω Sch x0 allows for a blowup near x 0 , all potentials in this class are integrable over compact sets, and so x 0 is still a regular point for the differential equation [9, footnote, p.67]; elsewhere the potentials are even more well-behaved.
1.1.The Solvability Complexity Index Hierarchy.Our exploration of the spectral computational problem -for both discrete and Schrödinger operators -continues a line of research initiated by Hansen in [20] and then further expanded in [5,4].This sequence of papers established the socalled Solvability Complexity Index (SCI) Hierarchy, which is a classification of the computational complexity of problems that cannot be computed in finite time, only approximated.Recent years have seen a flurry of activity in this direction.We point out [11,10] where some of the theory of spectral computations has been further developed; [31] where this has been applied to certain classes of unbounded operators; [2] where solutions of PDEs were considered; [6,7] where we considered resonance problems; and [12] where the authors give further examples of how to perform certain spectral computations with error bounds.
While precise definitions are provided below in Section 2, let us provide an informal overview.Computational problems can be classed as belonging to one (or more) of ∆ k , Σ k , Π k for k ∈ N as follows: ∆ k : For k ≥ 2, ∆ k is the class of problems that require at most k − 1 successive limits to solve.We also say that these problem have an SCI value of at most k − 1. Problems in ∆ 1 can be solved in one limit with known error bounds.Σ k : For all k ∈ N, Σ k ⊂ ∆ k+1 is the class of problems in ∆ k+1 that can be approximated from "below" with known error bounds.Π k : For all k ∈ N, Π k ⊂ ∆ k+1 is the class of problems in ∆ k+1 that can be approximated from "above" with known error bounds.By an approximation from "above" (resp."below") we mean that the output of the algorithm is a superset (resp.subset) of the object we are computing (this clearly requires that this object and its approximations belong to a certain topological space).It can also be shown that for k ∈ {1, 2, 3} we have In [4] the spectral computational problems for both B( 2 (N)) and for Schrödinger operators were addressed.Some of the results shown there include approximating σ(A) for A ∈ B( 2 (N)) We point out that showing that a problem belongs to Π 1 , Σ 1 or ∆ 1 is significant, as it shows that the computation can be done with certain guaranteed error bounds.
1.2.Periodic Operators.Schrödinger equations with periodic potentials have been the subject of study since the earliest days of quantum mechanics, perhaps most famously for the Bethe-Sommerfeld Conjecture [35], which states that the number of gaps in the essential spectrum is finite in dimensions ≥ 2. After more than seventy years this conjecture was finally proved in complete generality for Schrödinger operators in R d by Parnovski [26].In dimensions 2 and 3 the conjecture had already been proved by Popov and Skriganov [27] and Skriganov [34] respectively, and in dimension 4 by Helffer and Mohamed [21].
Beyond these results in mathematical analysis, since the 1990s interest in periodic problems has grown rapidly in the applied analysis and computational mathematics literature, partly driven by models of photonic crystals.These models are typically based on time-harmonic Maxwell equations or upon second order elliptic equations with periodic coefficients.Figotin and Kuchment [15,16,17] give particularly thorough analyses of some of these models, showing that already in these cases with piecewise constant coefficients the associated operators may possess an arbitrarily large number of spectral gaps.
For a periodic problem with some particular coefficients, often the first question of interest is whether it has any spectral gaps at all.Numerical methods may be used to obtain some preliminary evidence, and are almost always based on the Floquet-Bloch decomposition.The fact that the coefficients are often only piecewise continuous requires substantial effort to be given to adaptive meshing, see Giani and Graham [18], although the continuous variation of the quasi-momentum over the Brillouin zone means that some of the effort can be recycled from one quasi-momentum to the next.Despite the substantial computational cost, Floquet-Bloch techniques can sometimes be used to go beyond preliminary evidence, and have been combined with interval arithmetic to obtain algorithms which yield computer-assisted proofs of existence of spectral gaps for a wide class of problems with coefficients expressed in terms of elementary functions, see Hoang, Plum and Wieners [22].The literature for periodic problems which are not self-adjoint is much less extensive.For the ODE case, the work of Rofe-Beketov [30] is generally the starting point for any research on this topic.
The discrete case encompasses numerous different directions of research.Typical examples are Toeplitz and Laurent operators [19] or Jacobi operators [36].Another direction that has gained a great deal of attention in the last two decades is that of periodic discrete Schrödinger operators and generalizations thereof.A particular type known as almost Mathieu operator has been shown to exhibit rich spectral behavior (e.g. the spectrum can be a set of non-integer Hausdorff dimension, cf.[1,23]).Beyond this, the literature pertaining to operators that are either not tri-diagonal or not self-adjoint is limited.The best starting point would be the book of Trefethen and Embree [37].
Organization of the paper.In Section 2 we give a brief introduction to the main ideas of the SCI theory.Sections 3, 4 and 5 are dedicated to the proofs of Theorems 1.1, 1.2 and 1.3, respectively.Finally, in Section 6 we provide some numerical examples.

The Solvability Complexity Index Hierarchy
The Solvability Complexity Index (SCI) and the SCI Hierarchy provide a unified approach for understanding just how "difficult" it is to approximate infinite-dimensional problems (such as computing spectra) starting from finite-dimensional approximations.We start by setting the scene with a concrete example before providing precise abstract definitions.

Informal discussion & examples.
Consider the set Ω = B( 2 (N)) of all bounded operators on 2 (N).Let {e i } i∈N be the canonical basis.Then any element A ∈ Ω is represented by an infinite matrix.Denote by Λ the set of all entries in this matrix.Then one could ask: For any A ∈ Ω, is it possible to compute its spectrum σ(A) as the limit of a sequence of computations Γ n , where each Γ n has access to only finitely many elements of Λ and can only perform finitely many arithmetic computations?
Needless to say, the whole point here is that the algorithms Γ n are not tailored for this specific element A: they are meant to be able to handle any element A ∈ Ω.The convergence of Γ n (A) to σ(A) is made precise by realizing them as elements of the metric space M = (cl(C), d) which comprises all closed subsets of C endowed with an appropriate metric d (such as the Hausdorff metric).
In [20], Hansen showed that it is possible to compute σ(A) for any A ∈ Ω as above.However, rather than having algorithms Γ n with a single index n ∈ N, three indices were required, satisfying σ(A) = lim n3→∞ lim n2→∞ lim n1→∞ Γ n1,n2,n3 (A).The algorithms Γ n1,n2,n3 are given explicitly, and can be implemented numerically (though this raises a philosophical question about what it means to take successive limits numerically).In [4] it was proved that this is optimal: this computation cannot be performed with fewer than 3 limits, and hence we say that this problem has an SCI value of 3.
The SCI value strongly depends on Ω: intuitively, if Ω contains fewer elements, then devising a 'one-size-fits-all' algorithm should be easier.Indeed, if one considers Ω sa := {A ∈ Ω | A is selfadjoint} then the SCI value reduces to 2 and for Ω cpt := {A ∈ Ω | A is compact} it further reduces to 1.
The classification into SCI values can be further refined into a classification that takes into account error bounds.This is the so-called SCI Hierarchy which we describe below.2.2.Definitions.We formalize the foregoing example with precise definitions: Definition 2.1 (Computational problem).A computational problem is a quadruple (Ω, Λ, Ξ, M), where (i) Ω is a set, called the primary set, (ii) Λ is a set of complex-valued functions on Ω, called the evaluation set, (iii) M is a metric space, (iv) Ξ : Ω → M is a map, called the problem function.
Remark 2.2.In this paper it is often clear what M, Λ and Ξ are, and the important element is the primary set Ω.In this case we may abuse notation and refer to Ω alone as the computational problem.
Definition 2.3 (Arithmetic algorithm).Let (Ω, Λ, Ξ, M) be a computational problem.An arithmetic algorithm is a map Γ : Ω → M such that for each T ∈ Ω there exists a finite subset Λ Γ (T ) ⊂ Λ such that (i) the action of Γ on T depends only on {f (T )} f ∈ΛΓ(T ) , (ii) for every S ∈ Ω with f (T ) = f (S) for all f ∈ Λ Γ (T ) one has Λ Γ (S) = Λ Γ (T ), (iii) the action of Γ on T consists of performing only finitely many arithmetic operations on {f (T )} f ∈ΛΓ(T ) .
Definition 2.4 (Tower of arithmetic algorithms).Let (Ω, Λ, Ξ, M) be a computational problem.A tower of algorithms of height k for Ξ is a family Γ n1,n2,...,n k : Ω → M of arithmetic algorithms such that for all T ∈ Ω Ξ(T ) = lim Definition 2.5 (SCI).A computational problem (Ω, Λ, Ξ, M) is said to have a Solvability Complexity Index (SCI) of k if k is the smallest integer for which there exists a tower of algorithms of height k for Ξ.If a computational problem has solvability complexity index k, we write If there exists a family {Γ n } n∈N of arithmetic algorithms and N 1 ∈ N such that Ξ = Γ N1 then we define SCI(Ω, Λ, Ξ, M) = 0.
Definition 2.6 (The SCI Hierarchy).The SCI Hierarchy is a hierarchy {∆ k } k∈N0 of classes of computational problems (Ω, Λ, Ξ, M), where each ∆ k is defined as the collection of all computational problems satisfying: with the special class ∆ 1 defined as the class of all computational problems in ∆ 2 with a convergence rate: Hence we have that ∆ When the metric space M has certain ordering properties, one can define further classes that take into account convergence from below/above and associated error bounds.In order to not burden the reader with unnecessary definitions, we provide the definition that is relevant to the cases where M is the space of closed (and bounded) subsets of R d together with the Attouch-Wets (Hausdorff) distance (definitions of which can be found in Appendix A).These are the cases of relevance to us.A more comprehensive and abstract definition can be found in [4].

Banded Periodic Matrices
In this section we prove Theorem 1.1.We do this by defining an explicit algorithm and show that its output converges to the desired spectrum.The two computational problems we consider only differ in the primary set Ω, and are as follows: where {e i } i∈Z denotes the canonical basis of 2 (Z) and d H denotes the Hausdorff distance.We remind the reader that Ω per N,b is the class of operators on 2 (Z) whose canonical matrix representation has bandwidth b (i.e.A ij = 0 ∀|i − j| > b) and whose matrix entries repeat periodically along the diagonals with period N (here N, b ∈ N).Clearly, every A ∈ Ω per N,b defines a bounded operator on 2 (Z).Note that Ω per = N,b∈N Ω per N,b is the class of operators on 2 (Z) whose canonical matrix representation is banded and whose matrix entries repeat periodically along the diagonals.In the language of the Solvability Complexity Index, the two parts of Theorem 1.1 can be expressed as follows: • Part (i) amounts to proving that the computational problem for Ω per has an SCI value of 1 (or, equivalently, it belongs to ∆ 2 ).• Part (ii) amounts to showing that the computational problem for Ω per N,b belongs to ∆ 1 , i.e. it can be approximated with explicit error bounds; this is restated as Theorem 3.3 below.Remark 3.1.We note that the Hausdorff distance is only defined for non-empty sets, and it is finite only if the sets are bounded.Hence it is important to observe that for any A ∈ Ω per , the set σ(A) is both non-empty and bounded.Indeed, boundedness of the spectrum follows immediately from boundedness of A, while non-emptyness follows from the Floquet-Bloch theory described in Section 3.2.We discuss the metrics used in this paper in Appendix A.
Example 3.2.The class Ω per N,1 contains all Jacobi-type matrices of the form 3.1.Proof of Theorem 1.1(i).To prove Theorem 1.1(i) we assume to be known Theorem 1.1(ii).Theorem 1.1(ii) can be restated in the language of the SCI Hierarchy as follows: Theorem 3.3.For any fixed N, b ∈ N the computational problem for Ω per N,b can be solved in one limit with explicit error bounds, i.e.Ω per N,b , Λ, Ξ, M ∈ ∆ 1 .The proof of this theorem is contained in Subsection 3.4 below, after some preparatory work.First, however, we prove Theorem 1.1(i): Proof of Theorem 1.1(i).By Theorem 3.3, for every N ∈ N there exists a family of algorithms {Γ To clarify the meaning of d i we note that in Example 3.2 one would have d i = (. . ., 0, c i , a i , b i , 0, . . .).Loosely speaking, Pseudocode 1 first takes a finite section of A, searches it for periodic repetitions, and then defines a matrix B n ∈ Ω per N,b by periodic extension, to which Γ (N ) n can be applied.Because A is banded and periodic, this routine will eventually find its period: there exists n 0 ∈ N such that for all n > n 0 , N (as defined in the routine) is equal to the period of A and Finally, note that every line of Pseudocode 1 can be executed with finitely many algebraic operations on the matrix elements of A.
The following subsections are devoted to the proof of Theorem 3.3.The proof is constructive, i.e. we will provide an explicit algorithm that computes the spectrum of any given operator with explicit error bounds.

Floquet-Bloch transform.
Let N be as in the statement of Theorem 3.3.Given a vector We also introduce the symbol 2 per (N ) to denote the space of all N -periodic sequences (y k ) k∈Z , together with the norm y 2 Proof.This is standard, and the proof is omitted.

Transform and Properties of
where we have used the periodicity of A (see (3.4)) in the third line.The assertion now follows from the invertibility of U.
Remark 3.6.Observe that (A(θ)y) n = (A(θ)y) n+kN for any k ∈ Z so that A(θ)y ∈ 2 per (N ).Hence, A(θ) is an operator 2 per (N ) → 2 per (N ), which can be expressed as an N × N matrix.Note, however, that the numbers e −iθ n−j N A nj are not the matrix elements of A(θ) with respect to any basis.Indeed, Note that the sum in the last line is actually finite, because A is banded.Indeed, if the band width of A is less than the period, then the sum over j in (3.5) contains only one term.
Example 3.7.If A is a matrix with N = 1, i.e.A is a Laurent operator, then formula (3.5) yields a scalar function of θ given by Writing z := e iθ , we see that A(θ) = j∈Z z j A 0,j is given by the symbol of the Laurent operator.We thus recover the classical result that the spectrum of a Laurent operator is given by the image of the unit circle under its symbol (cf.[37,Th. 7.1]).
Next, we establish some elementary facts about the spectrum of a periodic operator.By standard results about the Floquet-Bloch transform, we have . Thus, an algorithm may be devised by determining the zeros of the map z → det(A(θ) − zI), θ ∈ [0, 2π].To this end, note that by definition of the determinant, one has where the coefficient functions p n (θ) are polynomials in the matrix entries A(θ) mn and hence analytic and periodic in θ.Hence they are bounded: Moreover, p N (θ) ≡ 1.
where we note that A ∞ = max{|A ij | | i, j ∈ Z} can be computed in finitely many steps. Proof.
From the mean value theorem it follows that for any differentiable function This follows from the Jacobi formula: for any square matrix M one has where cof(M ) denotes the cofactor matrix of M .Hence, Using Hadamard's inequality to bound the cofactor matrix, we obtain the bounds where the last two lines follow from the explicit formula (3.5).The bounds (3.6) and (3.7) imply and the claim follows.
3.4.Proof of Theorem 1.1(ii).We can finally prove Theorem 1.1(ii) which was restated equivalently as Theorem 3.3.First, we define the family of algorithms {Γ n } n∈N , where each of them maps Γ n : Ω per N,b → M (we recall that M is the space of all compact subsets of C endowed with the Hausdorff metric).It is easy to see that for any follows from Young's inequality), a quantity which can be computed in finitely many steps.Therefore, if we denote Remark 3.11.We emphasize that (3.8) can be computed in finitely many arithmetic operations on the matrix elements of A. Indeed, computing the radius R A consists of a finite sequence of multiplications and additions, as does the computation of each of the determinants det(A(θ Proof of Theorem 3.3 (equiv.Theorem 1.1(ii)).The proof has two steps.
Step 1: σ(A) is approximated from above by Γ n (A).For any set K ⊂ C we denote by This inequality implies that z n ∈ Γ n (A) as soon as 2n Note that the right-hand side of (3.9) is computable in finitely many arithmetic operations if N and b are known a priori.
Step 2: σ(A) is approximated from below by Γ n (A).Next we prove that Γ n (A) ⊂ B δ (σ(A)) for n large enough.We first note that, since det(zI − A(θ)) is a polynomial in z, it can be factored to take the form where z i (θ) are the zeros of z → det(zI − A(θ)) (note that for a characteristic polynomial the coefficient of the leading order term is always 1).From (3.10) we obtain the bound Let z n ∈ Γ n (A) be an arbitrary sequence.Then, by definition, Together, steps 1 and 2 imply that for any given δ > 0 one has both Γ n (A) ⊂ B δ (σ(A)) and Since the right-hand side of (3.12) is computable from N , b and the matrix elements of A in finitely many arithmetic operations, we conclude that the computational problem is in ∆ 1 .
Example 3.12.The algorithm from Definition 3.10 can easily be implemented in Matlab.An example calculation with bandwidth b = 1, N = 5 and A of the form (3.2) with yields the following output in the complex plane.Note that the output is a set of points on a discrete grid in the complex plane and hence the output looks 'fat'.As n is taken larger this set of points dwindles to just those points that lie in an ever decreasing neighborhood of the true spectrum.The Matlab implementation that produced Figure 2 is available online at https://github.com/frank-roesler/TriSpec.

Schrödinger Operators with Periodic Potentials
In this section we prove Theorem 1.2 regarding the spectrum of Schrödinger operators H = −∆+V with periodic potentials V : R d → C. Again, this is done by defining an explicit algorithm.We shall consider three computational problems which only differ in their primary set Ω, and are as follows (primary sets are defined immediately below): where −∆+V is meant to be defined on L 2 (R d ) with domain H 2 (R d ), and d AW denotes the Attouch-Wets metric, which is a generalization of the Hausdorff metric for the case of sets which may be unbounded (see Appendix A for a brief discussion).Note that the spectrum σ(−∆ + V ) is always non-empty in this case, so taking this metric makes sense.For p > d and M > 0, the primary sets are defined as follows: Note that by Morrey's inequality, every V ∈ Ω Sch p is continuous, and so the evaluation set Λ which comprises point evaluations of V , is well-defined.In the language of the Solvability Complexity Index, the two parts of Theorem 1.2 can be expressed as follows: • Part (i) amounts to proving that the computational problem for Ω Sch has an SCI value of 1 (or, equivalently, it belongs to ∆ 2 ).• Part (ii) amounts to showing that the computational problem for Ω Sch p,M belongs to Π 1 , i.e. it can be approximated from above with explicit error bounds.
The proof of Theorem 1.2 is contained in Subsection 4.5, and it relies on the following weaker theorem: Theorem 4.1.The computational problem for Ω Sch p can be solved in one limit: Note that this theorem is evidently weaker than Theorem 1.2(i) as the class of potentials Ω Sch p considered here is a strict subset of the class Ω Sch = p>d Ω Sch p considered in Theorem 1.2(i).The proof is constructive, i.e. we provide an explicit algorithm that computes the spectrum of any given operator with V ∈ Ω Sch p .Note that this problem is fundamentally different from the discrete problem (3.1)where we could directly access the matrix elements of the (discrete) operator.Instead, the evaluation set Λ gives access to the point values of the potential.Hence our task is to construct a sequence of algorithms {Γ n } n∈N , such that each Γ n computes its output from finitely many point evaluations of V using finitely many algebraic operations.
The proof of Theorem 4.1 is contained in Subsection 4.4.Prior to that, Subsection 4.1 is dedicated to the Floquet-Bloch transform, in Subsection 4.2 we approximate the potential V and in Subsection 4.3 we define a provisional algorithm.

4.1.
Floquet-Bloch transform.The Floquet-Bloch transform for Schrödinger operators with periodic potentials is well-studied.The following lemma is a collection of results in [28] (below, S(R d ) denotes the Schwartz space of rapidly decaying functions).
Then the following hold.
(i) U θ extends uniquely to a bounded operator on L 2 (R d ); where Moreover, the map θ → H(θ) is analytic and The spectral identity (4.2) follows from the unitarity of U and a straightforward calculation, noting that θ∈[0,2π] d σ(H(θ)) is closed by analyticity and periodicity in θ.

4.2.
Estimating and Approximating the Potential.The critical step is to compute approximations to the spectrum σ(−∆ + V ) using only finitely many pointwise evaluations of V .This is done using the Floquet-Bloch transform in conjunction with the Birman-Schwinger principle.The approximate potential is defined in (4.10) and the critical error bound is stated in Lemma 4.7.

4.2.1.
Birman-Schwinger principle.Let p > d, V ∈ Ω Sch p , θ ∈ [0, 2π] d and let H(θ) be the corresponding Floquet-Bloch operator as in (4.1).Expanding the operator square, H(θ) can be written as Let us choose the following decomposition of H(θ).We define Then clearly one has H(θ) = H 0 + B(θ).The auxiliary constant 1, which is added in H 0 and subtracted again in B(θ) was chosen for convenience, so that H 0 becomes a positive invertible operator.Note that B(θ) is relatively compact with respect to H 0 .For λ / ∈ σ(H 0 ) one has , where I denotes the identity operator on L 2 ((0, 1) 0 is invertible, and in that case This identity (sometimes called the Birman-Schwinger principle) implies that 4.2.2.Schatten class estimates.We now study the analytic operator valued function We choose a Fourier basis for L 2 ((0, 1) d ), that is, we choose a numbering N j → k j ∈ 2πZ d such that |k j | is monotonically increasing with j and set e j := e ikj •x .
We note that e j ∈ dom(H(θ)) for all θ and e j L 2 ((0,1) d ) = 1 for all j ∈ N. In this basis, the operators H 1 /2 0 , H 0 and θ • ∇ are all diagonal and one has H Therefore, we have Now the following lemma is easily proved. where −1 and • denotes the L 2 operator norm.
Proof.Let λ ∈ C \ σ(H 0 ) and note that C λ < +∞ by our choice of λ.Observe that simple geometric considerations lead to the bound Next we turn to the potential term in K(λ, θ), that is, the operator H This is easily treated by the ideal property of C s , since the operator |θ| 2 − 1 + V is bounded.Indeed, we have for every s > d where the last line follows from a similar calculation to (4.6) and the fact that H − 1 /2 0 = 1 (this follows from the matrix representation (4.4) and the fact that k 1 = 0).Lemma 4.4 (Lipschitz continuity of K).For every s > d and λ, µ where C λ is defined as in Lemma 4.3 and c V = 1 2 + πd Proof.First, note that Using (4.9) together with the resolvent identity for (λ − H 0 ) −1 one obtains Taking norms in C s and using the estimates from the proof of Lemma 4.3 (and in particular (4.7)) we obtain the bound Finally, applying Lemma 4.3 and using |θ + ϑ| ≤ 4πd This concludes the proof.
While Lemma 4.4 gives precise information about the dependence of the Lipschitz constant of K on all parameters, it will be useful for us to have a bound which is less precise but more explicit (and manifestly computable).
Proof.This follows immediately from Lemma 4.4 via rather crude estimates, noting that [8, Ch. 9.3]).4.2.3.Approximation of the potential.Next we study the matrix representation of the potential V in the Fourier basis {e j } j∈N and its approximations.First, we observe that in the Fourier basis {e j } j∈N one has e j , V e m = Vkm−kj where Vk denote the Fourier coefficients of V .Indeed, a direct calculation gives e j , V e m = (0,1) d e ikj •x V (x)e ikm•x dx = (0,1) d V (x)e i(km−kj )•x dx = Vkm−kj .Now, we want to build a computable, finite size approximation of the matrix (V jm ) = ( e j , V e m ).We start by approximating the Fourier coefficients Vk .Lemma 4.6.Let n ∈ N and define the lattice Proof.We write Then, comparing the two sums term by term, we have where Morrey's inequality was used in the third line (cf.(28) in the proof of [8, Th. 9.12]).Summing these inequalities over I n , we finally obtain Let us introduce the approximate Fourier coefficients for k ∈ 2πZ d and n ∈ N, Note that the V appr,n k can be computed in finitely many operations from the information provided in Λ. Lemma 4.6 applied to the function f (x) = V (x)e ik•x leads to the error estimate Next, we define the approximate potential matrix V appr,n := V appr,n km−kj m,j∈N .We remark that the approximated potential matrix cannot be computed in finitely many algebraic operations from the point values of V , because it has infinitely many entries.This issue will be addressed next.Lemma 4.7 (Main Error Bound).For N ∈ N let H N = Span{e 1 , . . ., e N } and let P N : L 2 ((0, 1) d ) → H N be the orthogonal projection.Moreover, define Then for every s > d one has where we recall that Proof.Again, we denote by • the L 2 (R d ) operator norm in this proof.We first treat the θ • ∇ term.By equations (4.4) -(4.5) we have and by (4.6) we can estimate the above as To estimate the next term, we denote W n := |θ| 2 − 1 + V appr,n and W := |θ| 2 − 1 + V for brevity, and note that the diagonal operators H − 1 2 0 and (λ − H 0 ) −1 commute with P N .Thus we have Let us first consider the second term on the right-hand side of (4.13).
where the third line follows from (4.11), with a similar calculation as in (4.8).To simplify notation, we collect all constants independent of λ, n, N into one and write Next we turn to the first term on the right-hand side of (4.13).We add and subtract P N W and use the triangle inequality to obtain Next we note that, by (4.4) -(4.5),H Finally, we employ (4.4) -(4.5) again to estimate the finite section error in (4.15).A straightforward calculation shows that . Using these bounds in (4.15), we finally obtain the error estimate Combining (4.12), (4.14) and (4.16) yields the assertion with constants

4.3.
A Provisional Algorithm.Armed with the error estimates from the last section, we are now able to define an algorithm based on the operator K(λ, θ).We first recall a result on Lipschitz continuity of perturbation determinants.
, † and define ) .(4.17) Note that for every N ∈ N, Γ z0 N (V ) can be computed from the information in Λ using finitely many algebraic operations (recall in particular the approximated Fourier coefficients (4.10)).Therefore, every Γ z0 N defines an arithmetic algorithm in the sense of Definition 2.3.Proposition 4.10 (Convergence of Provisional Algorithm).Let z 0 ∈ C and let V ∈ Ω Sch p .The following statements hold.
(i) For any sequence , G are explicit constants, which can be computed in finitely many operations from z 0 , δ, s, p, d and the a priori bound M for V W 1,p .
Proof.(i) Let z N ∈ Γ z0 N (V ) and assume that z N → z for some z ∈ C. We need to show that z ∈ σ(−∆+V ).Since z N ∈ Γ z0 N (V ), we have det p I −P N K appr n(N ) (z, θ i N )P N → 0 for some sequence {θ i N } N ∈N ⊂ [0, 2π] d .Then there exists a convergent subsequence (again denoted by θ i N ) converging to some some θ ∈ [0, 2π] d .We first note that due to Theorem 4.8 we have the N -independent determinant error bound ) is well-defined by Lemma 4.3 because p > d).We note that the exponential factor in the last line is uniformly bounded in N by some explicit constant c exp (cf.Lemmas 4.3 and 4.7).Using Lemma 4.7 with our choice n(N ) = N α and s = p, we obtain Note that each term on the right hand side has a negative power of N and therefore tends to 0 as −1 remains bounded as N → +∞ by our assumption that z / ∈ j (1 + |k j | 2 ).Therefore we have To conclude the proof of (i), if we note that continuity in (z, θ) and periodicity in θ imply Hence we have 1 ∈ σ(K(z, θ)) and thus z ∈ σ(−∆ + V ).
(ii) Conversely, denote ).We need to show that there exists a sequence Thus, using (4.20) and the main error estimate (4.19), we obtain To keep the notation simple, we collect all constants independent of N into a single constant G = G(p, d, V ) (note that |θ N | ≤ 2π √ d for all N ).This gives where the last line holds for N large enough (note that 1 2d − 1 2p > 0 by our assumption p > d and that C z N remains bounded as N → +∞ by our assumptions on z).Comparing (4.21) to (4.17), we see that N , the existence of which is guaranteed by the definition of Θ N .By Corollary 4.5 and Theorem 4.8 this implies where Turning to the finite approximation K appr n , this gives where G denotes the explicit and computable constant introduced in step (ii) above and we have used the bound which immediately implies the assertion, noting that |z − z N | < 1 N < ε as soon as N > ε −1 .4.4.Proof of Theorem 4.1.The provisional algorithm Γ z0 N defined in (4.17) is not sufficient to prove Theorem 4.1, because it does not approximate any eigenvalues lying in the set j (1 + |k j | 2 ).This is ultimately due to the decomposition (4.3) into H 0 and B(θ).In order to solve this issue, we define a second provisional algorithm based on a different decomposition.We define Then, obviously, H0 + B(θ) = H(θ) and by the Birman-Schwinger principle, Proof.The assertion is equivalent to the following equation having no solutions for any integers n i , ñi : Clearly the right-hand side of the above equation is always irrational or zero, while the left hand side is always nonzero rational.Thus no solution exists.Subsections 4.2 and 4.3 carry over trivially to the decomposition H0 , B(θ).Analogously to (4.17) we define the new algorithm Then we choose n(N ) := N α (with α as in Definition 4.9) and define where now From the proofs in Subsections 4.2 and 4.3 we immediately obtain the following convergence result.Proposition 4.13 (Convergence of Second Provisional Algorithm).Let z 0 ∈ C and let V ∈ Ω Sch p .The following statements hold.
(i) For any sequence as soon as N > max{ε −1 , N δ,z0 }, where N δ,z0 was defined in (4.18).Armed with the provisional algorithms (4.17) and (4.23), we are finally able to define the main algorithm: Combining Propositions 4.10 and 4.13, we can complete the proof of Theorem 4.1: Proof of Theorem 4.1.Let p > d.For any V ∈ Ω Sch p we need to show that Indeed, by Propositions 4.10 and 4.13, the following holds (a) For any sequence Attouch-Wets convergence now follows from a standard argument, cf.[ Proof of Theorem 1.2.We can finally prove Theorem 1.2.

Proof of part (i).
The a priori knowledge of p can be removed by slightly modifying the provisional algorithms (4.17) and (4.23).An algorithm that achieves SCI = 1 is obtained by modifying the determinant threshold in (4.17 Proof of part (ii).Fix R > 0. According to our choice of numbering {Z i } i∈N we have R , which is greater than 0 by Lemma 4.11.This choice implies that removing δ R -neighborhoods of the free spectra does not actually restrict our domain of computation: for all R > 0. Now let ε > 0. By Propositions 4.10(iii) and 4.13(iii) we have We recall that • the lower bound 4R 2 ensures that B R (0) is covered by the search regions Q Zi ; • the lower bound ε −1 ensures that for every z ∈ σ(−∆ + V ) there exists To conclude the proof of part (ii) (i.e., that Ω Sch p,M ∈ Π 1 ), we need a set X k as in Definition 2.7.To this end, for V ∈ Ω Sch p,M and k ∈ N, choose N (k) := max{4k 2 , 2 k , N k } and define Then by the definition of the Attouch-Wets distance d AW (cf.Definition A.2) one has where the third line follows from the definition of X k .Moreover, by (4.24) (with ε = 2 −k and R = k) we have for all k ∈ N.These facts, together with step (i) imply Ω Sch p,M ∈ Π 1 .
Remark 4.15.In view of Theorems 3.3 and 4.1 it is natural to ask whether one might have Ω Sch p,M ∈ Σ 1 (and therefore Ω Sch p,M ∈ ∆ 1 ).This is indeed a nontrivial open problem.Recalling the proof of Theorem 3.3 (in particular (3.11)), the error bound for the approximation "from below" used the fact that det(zI − A(θ)) is a polynomial.In the proof of Theorem 4.1 on the other hand, the function det p (I − K(z, θ)), whose zeros must be approximated, is only known to be analytic.Obtaining Σ 1 classification amounts to obtaining explicit upper bounds on the width of the zeros of this analytic function.These can not be deduced in any straightforward way from the values of the potential V .

Sometimes Two Limits Are Necessary
In Section 2 we described a complicated construction called towers of algorithms, where more than one successive limit is required to correctly perform certain computations.In this section we exhibit this phenomenon first hand: we prove Theorem 1.3 which is rephrased in the language of SCI as Theorem 5.1 below.This theorem shows that there exists a class of potentials (which are less smooth, yet can be evaluated at any point) for which there do not exist algorithms that can approximate the associated spectral problem in a single limit.That is, SCI > 1.We are able, though, to construct an arithmetic algorithm which converges by taking two successive limits.That is, SCI = 2.This class contains potentials which are allowed to have a singularity at a single point x 0 ∈ (0, 1), but are otherwise smooth: We emphasize that Ω Sch x0 contains only real-valued functions and that the pointwise evaluation V → V (x) is well-defined for all x ∈ R.Moreover, we remark that by [28,Th. XIII.96] the Schrödinger operator −∆ + V is well-defined and selfadjoint with domain H 2 (R).We shall therefore consider the computational problem and shall prove Theorem 5.1.The computational problem (5.1) has SCI = 2 (equivalently, it belongs to ∆ 3 \ ∆ 2 ).
The proof of Theorem 5.1 has two parts.To prove SCI ≤ 2 we construct an explicit tower of algorithms that computes the spectrum in two limits (cf.Definition 2.4).The proof of SCI > 1 is by contradiction.We assume the existence of a sequence {Γ n } n∈N with Γ n (V ) → σ(−∂ 2 x + V ) for every V ∈ Ω Sch x0 and via a diagonal process construct a potential V ∈ Ω Sch x0 such that Γ n (V ) σ(−∂ 2 x + V ), yielding the desired contradiction.5.1.Lemmas.We first collect some technical lemmas that will be necessary for the proof.Lemma 5.2.Denote by W (H) the numerical range of an operator H. Let V ∈ Ω Sch x0 .There exists where the last line follows from periodicity of V .Moreover, we have where the second line follows from the Sobolev embedding W 1,1 → L ∞ (cf.[8,Th. 8.8]) and the third follows by the product rule and Hölder's inequality.Summing (5.2) over i we obtain where C = 2C sup i∈Z χ i W 1,∞ is finite by choice of the functions χ i .Turning to the numerical range, we have for φ ∈ C ∞ 0 (R) , where the second line follows from (5.3).Finally, let δ := (4C and the proof is complete.
Proof.Define the sequence of test functions We note that φ n L 2 (R) → 1 as n → +∞.Then if Ṽ ∈ Ω Sch x0 with Ṽ ≤ 0 we have The assertion follows by letting n → +∞.
The next lemma is needed to construct the tower of algorithms that will compute the spectrum of elements in Ω Sch x0 .
Lemma 5.4.Let V ∈ Ω Sch x0 and for n ∈ N let ρ n ∈ W 1,∞ (R) be a periodic function such that ) where z / ∈ R and the right hand side is defined whenever (ρ n − 1)V (−∂ 2 x + V − z) −1 L 2 →L 2 < 1.We will prove that the right hand side of (5.5) converges to 0 in the norm resolvent sense.To this end let u ∈ L 2 (R) and compute where Hölder's inequality was used in the second line, periodicity of (ρ n − 1)V was used in the third line and the Sobolev embedding H 1 (R) → L ∞ (R) was used in the fouth line.Combining eqs.(5.5) and (5.6) we find that , which immediately implies norm resolvent convergence.Finally, an application of [29, Thms.VIII.[23][24] yields the desired spectral convergence.
Step 1: construction of an adversarial potential V .Assume for contradiction that there exists a sequence of algorithms {Γ n } n∈N such that for any We now describe a process that defines an "adversarial" potential V for which the sequence {Γ n } n∈N will necessarily fail.We begin with an example construction that illustrates how a potential can be obtained that "fools" Γ n for a single n.This construction is then iterated below to obtain a potential whose spectrum is not approximated by the entire sequence {Γ n } n∈N .For the sake of definiteness we make the arbitrary choice x 0 := 1 4 , though the proof does not depend on this choice.Moreover, we note that by Hölder's inequality one has f With the notation from Lemma 5.3, let V 1 be the periodic square well potential Our assumptions about Γ n imply that there exists n 1 ∈ N such that inf Re(Γ n (V )) < − 1 2 for all n ≥ n 1 .By the definition of an algorithm (Definition 2.3), Γ n1 (V ) depends only on finitely many elements of Λ, i.e. finitely many point values V (x i ), i ∈ {1, . . ., m n1 }.We may assume without loss of generality that the sets {x 1 , . . ., x mn } are growing with n.
For later reference, let ρ be a smooth, compactly supported function on R such that Next, define a new potential Ṽ1 by "thinning out" V 1 around the points x i .More concretely, let δ > 0 (to be determined later) and define l 1 := min |x i − x j | i, j ∈ {0, . . ., m n1 }, x i = x j (note that the point x 0 appears in the definition of l 1 ) and let Then by construction we have Ṽ1 and thus d AW Γ n1 ( Ṽ1 ), σ(−∂ 2 x + Ṽ1 ) ≥ 1 4 .We remark that Ṽ1 ≡ 0 in a neighbourhood of x 0 .The constructions outlined above define an iterative process that yields a sequence of smooth potentials {V k } k∈N and { Ṽk } k∈N .We outline the details below.Fix η > 0 to be determined later and initialize Ṽ0 ≡ 0.
Step 2: Construction of a tower of algorithms.We conclude the proof of Theorem 5. where all limits are taken in Attouch-Wets distance.This completes the proof.
Note that the two limits in Step 2 above cannot be swapped, because it is unclear how Γ m (ρ n V ) behaves when ρ n V converges to a non-smooth function.

Numerical Results
To illustrate our abstract results, we implemented a version of algorithm (4.14) in one dimension in Matlab.In this section we show the results of this implementation and compare them against known abstract and numerical results.
In order to obtain an implementation with adequate performance, we fixed a box Q in the complex plane and then computed the quantity where the numbers N , n, C were treated as independent parameters.Moreover, the spectral shift, which was chosen to be 1 in (4.3) and 2 in (4.22) can be fixed to be any point z 0 outside the box Q.
The routine is illustrated by the following pseudocode.where µ ∈ C is a constant and λ denotes the spectral parameter.This equation was first studied in [24] in the context of vibrating membranes and has been studied extensively since (see [25,Ch. 5] for a discussion).Figure 4 shows the output of our implementation for various values of µ.In the case of a real-valued potential (top panel of Figure 4) our algorithm produces the expected band-gap structure, with one gap showing around λ = 10 and another around λ = 40.In the case of purely imaginary µ, the theory of PT -symmetric operators can be used to prove abstract results about the possible shape of the spectrum [32].A comparison between the middle panel of Figure 4 and [32,Fig. 2] shows agreement between the theoretical results and the output of our algorithm.Finally, the bottom panel in Figure 4 shows the output when µ has both a real and an imaginary part and the PT symmetry is broken.
Finally, we note that our method is naturally immune to the common problem of spectral pollution.By definition, spectral pollution occurs when there exist sequences z n ∈ Γ n (V ) such that {z n } n∈N has an accumulation point outside σ(−∆ + V ).This effect appears in the approximation of Mathieu's equation if a naive finite difference scheme is used on a truncated domain (see Figure 6).For  arbitrarily fine discretization and arbitrarily large truncated domain, the method is sensitive to small perturbations of the domain and can yield a set which is far from the correct spectrum.

Definition 2 . 7 (
The SCI Hierarchy (Attouch-Wets/Hausdorff metric)).Consider the setup in Definition 2.6 assuming further that M = (cl(R d ), d) where d = d AW or d = d H . Then for k ∈ N we can define the following subsets of ∆ k+1 :

k=0
|y k | 2 .Note that 2 per (N ) is canonically isomorphic to the Euclidean space R N .The following lemma is easily proved by direct computation.Lemma 3.4 (Properties of U θ ).The map U θ defined in (3.3) has the following properties.

2 per 1 = 2 per
(N ) is finite-dimensional, while e −iθ n−j N A nj (n, j ∈ Z) are infinitely many numbers.As noted earlier, we have 2 per (N ) ∼ = R N .The identification can be made via the basis e per .(e per n ) j = δ (j mod N ), n for n ∈ {0, . . ., N − 1} (Kronecker symbol).In this basis, the matrix elements of A(θ) become A(θ) mn = e per m , A(θ)e per n

Figure 2 .
Figure 2. Output of Γn(A) in the complex plane for different values of n.

Lemma 4 . 3 (
Schatten bound for K).For every s > d the operator K(λ, θ) belongs to the Schatten class C s and one has

1 2d − 1 2p)
) and(4.23).If the cutoff N −( is replaced by log(N ) −1 and the discretization width n(N ) = N α is replaced by n(N ) := e N , the proof of Propositions 4.10 and 4.13 is valid independently of the value of p (cf. eq.(4.21)).

Figure 3 .
Figure 3. Sketch the first few iterations V k , Ṽk .The red dots represent the points xi, which represent the information in ΛΓ n .

Figure 6 .
Figure 6.Spectral approximation from finite difference scheme on truncated domain [−100.7,100.7].The spurious eigenvalue in the gap for µ = 10 can appear for arbitrarily fine discretization.
∈ Ω per and define a new family {Γ n } n∈N by the following pseudocode.Definition of {Γ n } n∈N on Ω per for n ∈ N do For i ∈ {−n, . . ., n} define d i := (a i,i−n , a i,i−n+1 , . . ., a i,i , . . ., a i,i+n−1 , a i,i+n ) if ∃p < n s.t.d i+p = d i ∀i ∈ {p, . . ., n − p} then N