Computing scattering resonances

. The question of whether it is possible to compute scattering resonances of Schrödinger operators – independently of the particular potential – is addressed. A positive answer is given, with the potential merely required to be C 1 and have compact support. The proof is constructive, providing a universal algorithm which only needs to access the values of the potential at any requested point. Numerical examples are provided and compared with known results


Introduction and Main Result
This paper provides an affirmative answer to the following question:

Does there exist a universal algorithm for computing the resonances of Schrödinger operators with complex potentials?
To the authors' best knowledge this is the first time this question is addressed.Furthermore, the proof of existence provides an actual algorithm (that is, the proof is constructive).We test this algorithm on some standard examples, and compare to known results.
The framework required for this analysis is furnished by the Solvability Complexity Index (SCI), which is an abstract theory for the classification of the computational complexity of problems that are infinite-dimensional.This framework has been developed over the last decade by Hansen and collaborators (cf.[14,4,5]) and draws inspiration from the seminal result [10] on solving quintic equations via a tower of algorithms.We therefore emphasize that ours is an abstract result in analysis, not in numerical analysis.
1.1.Quantum Resonances.Let us first define what a quantum resonance is.Let q : R d → C be compactly supported, let H q := −∆ + q be the associated Schrödinger operator in L 2 (R d ) and let χ : R d → R be some compactly supported function with χ ≡ 1 on supp(q).It follows from the explicit form of the free fundamental solution (cf.eq.(1.1) below) that the map is an analytic operator-valued function on C \ {0}.We define: Definition 1.1 (Resonance).A resonance of H q is defined to be a pole of the meromorphic operatorvalued function z → (I + q(−∆ − z 2 ) −1 χ) −1 .
Definition 1.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 SCI(Ω, Λ, Ξ, M) = k.
In the present article, our computational problem is made up of the following elements (to be specified more precisely in Section 1.3): (i) Ω is a class of Schrödinger operators H q with potentials q which have a common compact support and a uniform bound in C 1 , (ii) Λ is the set of all pointwise evaluations of q, as well as pointwise evaluations of the Green's function associated with the Helmholtz operator −∆ − z 2 , (iii) M is the space cl(C) of all closed subsets of C equipped with the Attouch-Wets metric (which is a generalization of the Hausdorff metric to the case of unbounded sets), (iv) Ξ : Ω → M is the map that associates to a particular Schrödinger operator its set of resonances, and we denote it by Res(H q ).
We show that for this computational problem there exists a tower of height 1, i.e. there exists a family of algorithms {Γ n } n∈N such that Γ n (H q ) → Res(H q ) as n → +∞ for any H q ∈ Ω, where the convergence is in the sense of the Attouch-Wets metric [3], generated by the following distance function: Definition 1.6 (Attouch-Wets distance).Let A, B be closed sets in C. The Attouch-Wets distance between them is defined as Remark 1.7.It can be shown (cf.[20,Prop. 2.8]) that a sequence of sets A n ⊂ C converges to A in Attouch-Wets metric, if the following two conditions are satisfied where H (1) ν denotes the Hankel function of the first kind.For Im(z) > 0, G(x, z) is the fundamental solution to the free Helmholtz operator −∆ − z 2 (cf.[21,Ch. 22]).We define the evaluation set Λ to be Then the quadruple (Ω M , Λ, Res(•), cl(C)) poses a computational problem in the sense of Definition 1.2.The main result of the present article is the following.
Theorem 1.8.Resonances can be computed in one limit: We prove this theorem by explicitly constructing an algorithm which computes the set of resonances in one limit.This algorithm can be implemented numerically; some numerical experiments are provided in Section 5.
Remark 1.9.(i) We note that there exist examples of computational spectral problems for which SCI ≥ 2, even in the selfadjoint case (cf.[4,Th. 6.5]).Whenever SCI ≥ 2 it is impossible to have error control (this can be shown using a straightforward diagonal-type argument).
Identifying situations in which SCI = 1 is therefore of particular interest.(ii) If we drop the restriction supp(q) ⋐ Q M from the definition of Ω, i.e. if the size of supp(q) is not assumed to be known a priori, then Theorem 1.8 implies SCI ≤ 2, but it is no longer clear that SCI = 1.Indeed, let us provide a sketch of the proof: Sketch of proof: Choose a sequence 0 ≤ M k → +∞ and for each fixed k let {Γ k,n } n∈N be a sequence of algorithms as in Theorem 1.8, which computes the resonances of Ω M k in one limit.Given a compactly supported potential function q, whose support is not known, choose a smooth cutoff function Then from Theorem 1.8 we know that {Γ k,n } n∈N computes the set of resonances of H ρ k q in one limit.As soon as supp(q) ⋐ Q M k −1 , one has ρ k q ≡ q and the sequence {lim n→∞ Γ nk (H ρ k q )} k∈N will be constant in k, and hence convergent.This proves that lim k→∞ lim n→∞ Γ k,n (q) = Res(H q ), that is, the set of resonances can be computed in two limits.
The proof of Theorem 1.8 is divided into several steps.First, we obtain quantitative resolvent norm estimates for the operator K(z) := q(H q − z 2 ) −1 χ from Definition 1.1.These are then used to bound the error between K(z) itself and a discretized version K n (z), obtained by replacing the potential q by a piecewise constant approximation.Finally, the poles of (I + K(z)) −1 are identified through a thresholding of the discretized operator function Organization of the paper.Section 2 contains a short discussion of Definition 1.1 and meromorphic continuation.In Section 3 we prove some estimates for convergence of finite-dimensional approximations of linear operators, which are then used in Section 4 to construct an explicit algorithm which computes resonances in one limit, thereby proving Theorem 1.8.Section 5 is dedicated to numerical experiments.In Appendix A we review some properties of the Green's function G(x, z) introduced in (1.1).

Analytic Continuation
We use this section for a more detailed discussion of Definition 1.1 and to fix some notations and conventions.For the sake of self containedness, we prove the existence of z → (I +q(−∆−z 2 ) −1 χ) −1 as a meromorphic operator-valued function on the domain This result follows from the classical Analytic Fredholm Theorem (cf.e.g.[18, Sec.VI.5]) Theorem 2.1 (Analytic Fredholm Theorem).Let D ⊂ C be open and connected and let F : D → L(H) be an analytic operator-valued function such that F (z) is compact for all z ∈ D.Then, either , the residues at the poles are finite rank operators, and if z ∈ S then ker(I + F (z)) = {0}.
Next, recall that Q M denotes the cube of edge length M in R d centered at the origin.Let χ := χ QM be the indicator function of Q M .Note that the operator-valued function z → q(−∆ − z 2 ) −1 χ is an analytic function on C ext \ {0}.This follows from the explicit representation of the free fundamental solution (1.1) (cf.Remark A.2).

Lemma 2.2. The function
Moreover, the residues at the poles are finite rank operators.
Proof.The operator q(−∆ − z 2 ) −1 χ is compact by the Fréchet-Kolmogorov theorem and the inverse The above observations lead us to study the spectrum of the compact operator Since the integral kernel for the free resolvent is given explicitly by (1.1) as an analytic function of z ∈ C ext \ {0}, we have an explicit representation of (2.1) as an integral operator on L 2 (R d ):

Abstract Error Estimates
We recall that the resonances of H q = −∆ + q are defined to be the poles of In this section we prove general, abstract, estimates for approximations of families of linear operators.These are largely independent of the rest of this paper and will be applied in the proof of Theorem 1.8.Abusing notation, our generic abstract analytic operator family is denoted K(z).
Let H be a separable Hilbert space and denote by L(H) the space of bounded operators on H. Let H n ⊂ H be a finite-dimensional subspace, P n : H → H n the orthogonal projection and K : C ext → L(H) continuous in operator norm.Moreover, let K n : C ext → L(H n ) be analytic for every n ∈ N. Assume that for any compact subset B ⊂ C ext there exist a sequence a n ↓ 0 and a constant C > 0 such that , where we use the convention that Proof.Whenever the left hand side is non-positive the assertion is trivially true, so we may assume w.l.o.g. that 1 In this case, the assertion follows by a Neumann series argument, as follows.We have , the second factor in (3.4) is invertible by the Neumann series and Hence, Ca n for any n ∈ N. It remains to replace the L(H) norm on the left hand side by the L(H n ) norm.This follows from Claim 3.2.This completes the proof.Claim 3.2.We have for all z for which both operators are boundedly invertible.
Proof.If −1 ∈ σ(K(z)), then unless −1 ∈ σ(P n K(z)P n ), we have We now argue by contradiction.If we had (I + P n K(z)P n ) −1 L(H) < 1 2Can , then we would have (I + P n K(z)P n ) −1 (K(z) − P n K(z)P n ) L(H) < 1 and I + K(z) would be invertible by the Neumann series contradicting our assumption that −1 ∈ σ(K(z)).Thus we must have (I + . Now let us turn to the case where −1 / ∈ σ(K(z)) and ( Can .The same calculation as in the proof of Lemma 3.1 shows that from which it follows easily that 1 2Can ≤ (I + P n K(z)P n ) −1 L(H) .3.2.An Abstract Algorithm For Computing Poles.We now demonstrate how the the assumptions (3.1)-(3.3)allow us to construct an abstract algorithm that computes the poles of I +K(z) −1 .
By an abstract algorithm we mean a sequence of subsets of C ext , which is constructed from K n and which converges in Attouch-Wets metric to {z ∈ C ext | − 1 ∈ σ(K(z))}.Note that this is not yet an arithmetic algorithm in the sense of Definition 1.3, since the sets are not computed from a finite amount of information in finitely many steps.Let B ⊂ C ext be compact and define the exponentially fine lattice L n := e − 1 an (Z + ıZ) ∩ B. Since we assume that a n is explicitly known and K n (z) can be computed in finitely many steps, we can define the set Moreover, note that by [4, Prop.10.1], determining whether ( an can be done with finitely many arithmetic operations on the matrix elements of K n (z) for each z ∈ L n .
an and hence by Lemma 3.1 (with the convention that ( It follows that (I + K(z n )) −1 L(H) → +∞ as n → +∞ and hence I + K(z 0 ) is not invertible (this follows by yet another Neumann series argument, together with norm continuity of K).Hence z 0 is a pole.
II. Spectral inclusion.Assume now that z is a pole, i.e. −1 ∈ σ(K(z)).Our reasoning will have the structure with a quantitative estimate in each step.To this end, note first that if −1 ∈ σ(K(z)) for some z ∈ B, then there exist ν, c, ε > 0 (independent of n) such that for all ζ in a ε-neighborhood of z, Indeed, since all singularities of (I + K(z)) −1 are of finite order by the analytic Fredholm theorem, this follows from the Laurent expansion of meromorphic operator valued functions.
It follows from (3.5) that for any z n such that |z − z n | ≤ e − 1 an one will have for all large enough n.We conclude that for any pole z there exists a sequence z n ∈ L n such that z n → z as n → +∞ and ( Can for all but finitely many n ∈ N. Next, Lemma 3.3 shows that (I + P n K(z n )P n ) −1 L(H) > 1 2Can .Studying this norm further, we have and thus Hence, as soon as a n < 1 2C , we have ( . We conclude that if z is a pole, then there exists z n ∈ L n such that (n large enough).A similar reasoning as in Lemma 3.1 (using (3.2)) shows that now and rearranging terms, together with (3.6), gives and therefore z n ∈ Θ B n (K) for large enough n.The assertion about Attouch-Wets convergence now follows from Remark 1.7.

Definition of the Algorithm
4.1.Error Estimates.In this section, we will apply the abstract results of Section 3 to our resonance problem.To this end, define K(z) = q(−∆ − z 2 ) −1 χ.We write K for the integral kernel of K(z) to simplify notation.Recall that K was given by (2.2) and supp(K) ⊂ Q M × Q M .We will construct an operator approximation K n of K, which satisfies (3.1)-(3.3)and in addition (H1) The matrix elements of K n can be computed in finitely many steps from a finite subset Λ n ⊂ Λ (cf.eq. ( 1.2) and Def.1.3); (H2) The convergence rate a n is explicitly known (i.e. the sequence a n can be used to define the algorithm).
To this end, let us define H n , P n as follows.
Furthermore, we have to make a concrete choice for the approximation K n .An obvious choice is the integral kernel i.e. a piecewise constant approximation of K(•, •) which can be computed from the values of K on the lattice n −1 Z d (in dimensions larger than one, the fundamental solution G has a singularity at x = y.Hence, we put K n := 0 for i = j in this case).
We will now show that the operators K, K n satisfy eqs.(3.1)-(3.3).To streamline the presentation, we will restrict ourselves to d ≥ 3 in our computations, the cases d ≤ 2 being entirely analogous with minor changes in the formulas.Constants independent of n will be denoted C and their value may change from line to line.Proof of (3.3).Using the definitions (4.1)-(4.2),we have where Note that P y n P x n K(x, y) simply yields a step function approximation of K(x, y) like (4.2), but in dimension 2d.We conclude by applying Young's inequality [23,Th. 0.3.1], that where Thus, all we have to do is estimate the L ∞ -L 1 difference between K and its projection onto step functions.To this end, fix x ∈ Q M , let ε > 2 n and decompose the integrals as follows The integral over B ε (x) can be estimated by Bε(x) 2|K(x, y)| dy, while for the remaining integral we can use the fact that the derivative of K is bounded, as follows.Let j ∈ 1 n Z d be such that x ∈ S n,j , see Figure 1 Figure 1.Sketch of the geometry in the calculation leading to (4.5).The sum over i includes all cells whose nodes are outside the dashed ball centered at j.
Summing over i, we finally obtain (cf. Figure 1) where the fifth line follows from (A.2), in the appendix, and the bound q C 1 ≤ +∞.Using (4.5) in (4.4), we conclude that where in the last line we have used (A.1) and the boundedness of q again.With an analogous calculation for sup y∈R d R d |K(x, y) − P y n P x n K(x, y)| dx (which we omit here), and recalling that η n was defined by (4.3), we conclude that for all ε > 0 Choosing ε := n − 1 d+1 , we conclude that Remark 4.1.Note that the constants C, C ′ all depend on the spectral parameter z, but are bounded for z in compact subsets of C ext , because K depends continuously on z.
Proof of (3.2) and (H1).An orthonormal basis of H n is given by the functions f, e j L 2 e j in this basis.It is then easily seen that in this basis K n has the matrix elements Note that this proves (H1): The matrix elements of K n can be calculated in finitely many arithmetic operations from the finite set Λ Similarly, it can be seen that the matrix elements of P n K| Hn in this basis are given by where we have introduced the notation • ij for the mean value on S n,i ×S n,j .Let f = j f j e j ∈ H n .From the above, and Young's inequality, we conclude that where ηn := max sup Hence, we have reduced the problem to estimating these ℓ ∞ -ℓ 1 differences.This can be done similarly to (4.4), by separating ) into an ε-region around i = j and the rest: where we have used (A.2) and the C 1 -boundedness of q in the last line.To estimate the last term on the right hand side, note that where ω d denotes the volume of the unit sphere in R d .Note that the above calculation is uniform in i, because q is bounded.Plugging (4.8) into (4.7),we arrive at Finally, swapping i and j will give an analogous estimate and we can conclude that ηn → 0 with rate a n = n − The explicit rates obtained in (4.6) and (4.9) prove that our approximation scheme satisfies (H2).
4.2.The Algorithm.It remains to extend the algorithm Θ B n from a single compact set B ⊂ C ext to the entire complex plane.This is done via a diagonal-type argument.Next, we define our algorithm as follows.We let for fixed k and since the {B (k) j } form a tiling of C ext , it follows that Γ n (q) → Res(q) in Attouch-Wets metric.The proof of Theorem 1.8 is complete.

Numerical Results
Software to compute resonances has been in existence for decades [19,9,1].The authors of [7] recently proposed a collection of MATLAB codes to compute resonance poles and scattering of plane waves efficiently ("MatScat", cf.[6]).In this section we compare the results of our algorithm to that of MatScat.
In order to study the actual numerical performance of our algorithm, we coded a MATLAB routine for the one-dimensional case with supp(q) ⊂ [a, b] (for some known a < b), which computes the set where the region B in the complex plane, the lattice distance of L n and cutoff threshold C were treated as independent parameters.Comparison of results.Figures 4 and 5 show the output of MatScat (black dots) versus the output of our algorithm (blue regions) for a Gaussian well and trapping potential, respectively.As the plots show, there is agreement between the two.Limitations.As mentioned before, MatScat has been developed with the goal to create an efficient algorithm to compute resonances fast.Indeed, the computation of the black dots in Figure 4 takes less than a second, while computing the regions with our algorithm takes several hours on a personal computer.We stress that our MATLAB code was written mainly for illustration purposes and that

1. 3 .
Main Result.Let d ∈ N, fix M > 0 and let Ω M denote the class of Schrödinger operators by the Neumann series.Hence, the claim follows from the analytic Fredholm theorem, together with Remark A.2 in the appendix.

Figure 2 .
Figure 2. Tiling of the complex plane

Figure 4 .
Figure 4. Comparison of the result of [6] (black) and our algorithm (blue) for a Gaussian well supported between −1 and 1.The chosen parameter values are: n = 100; threshold for resolvent norm: C = 200; number of lattice points in the shown region of the complex plane: M × 4M = 1000 × 4000.

Figure 6 .
Figure 6.Numerical artefacts for large real part of z.Top: Output of our algorithm for Gaussian well potential on the interval [−1, 1] with n = 15.Bottom: Output for the same problem with n = 30.The locations of the spurious peaks agree with the bound | Re(z)| ∼ πn in each case.
Remark 4.2.Note again that the constants C, C ′ depend on z, but are bounded for z in compact subsets of C ext , since K depends continuously on z.Proof of (3.1) and (H2).Estimate (3.1) in fact follows from (3.3) and (3.2).Indeed, writing K n and K as block operator matrices w.r.t. the decomposition H = H n ⊕ H ⊥ n , we have