Guaranteed resonance enclosures and exclosures for atoms and molecules

In this paper, we confirm, with absolute certainty, a conjecture on a certain oscillatory behaviour of higher auto-ionizing resonances of atoms and molecules beyond a threshold. These results not only definitely settle a more than 30 year old controversy in Rittby et al. (1981 Phys. Rev. A 24, 1636–1639 (doi:10.1103/PhysRevA.24.1636)) and Korsch et al. (1982 Phys. Rev. A 26, 1802–1803 (doi:10.1103/PhysRevA.26.1802)), but also provide new and reliable information on the threshold. Our interval-arithmetic-based method allows one, for the first time, to enclose and to exclude resonances with guaranteed certainty. The efficiency of our approach is demonstrated by the fact that we are able to show that the approximations in Rittby et al. (1981 Phys. Rev. A 24, 1636–1639 (doi:10.1103/PhysRevA.24.1636)) do lie near true resonances, whereas the approximations of higher resonances in Korsch et al. (1982 Phys. Rev. A 26, 1802–1803 (doi:10.1103/PhysRevA.26.1802)) do not, and further that there exist two new pairs of resonances as suggested in Abramov et al. (2001 J. Phys. A 34, 57–72 (doi:10.1088/0305-4470/34/1/304)).

In this paper, we confirm, with absolute certainty, a conjecture on a certain oscillatory behaviour of higher auto-ionizing resonances of atoms and molecules beyond a threshold. These results not only definitely settle a more than 30 year old controversy in Rittby et al. (1981 Phys. Rev

Introduction
Reliable and precise information on the location of resonances is very hard to obtain. While numerical approximations are widely used in physics, so far there has been no way to show that they produce results near, or not near, true resonances. The reason is that computations of complex eigenvalues in the presence of continuous spectrum are not backed up by any convergence results. This paper presents a new method that, for the first time, permits one to locate resonances with absolute certainty and high accuracy and, at the same time, to show that numerical approximations fail to lie near true resonances. We provide new and reliable information on the oscillatory behaviour of the real parts of certain resonance strings and on the threshold beyond which it occurs.
The key ingredient in our method is interval arithmetic. It allows us to carry out every computational step with absolute accuracy by operating on intervals rather than on numbers. Remarkably, this theoretical idea has had convincing impact in different practical physical applications recently: to control the stability of difficult nonlinear systems in robotics to navigate a sailboat autonomously over a distance of 100 km (see [1]); to perform rigorous global optimization of impulsive planet-to-planet transfer (see [2]) or to rigorously govern the long-term stability in particle accelerators (see [3]).
In this paper, we demonstrate the efficacy of interval approaches for the computation of resonance enclosures and exclosures with absolute certainty. The power of our method is substantiated by the fact that it can be applied to definitely settle a more than 30 year old controversy in [4,5] which could not be resolved by any other method before.
In connection with auto-ionizing resonances of atoms and molecules lying above a ionization threshold, Moiseyev et al. [6] studied resonances of the Sturm-Liouville problem (1.1) A first resonance was suggested to lie near 2.124 − 0.0185i (with μ/h 2 set to 1); moreover, one bound state was proposed to lie near 0.5. The resonance was found by complex scaling of the self-adjoint Hamiltonian and approximation using a variational principle with 10 real Gauss-type basis functions for the scaled Hamiltonian. Because the latter is no longer self-adjoint, the authors pointed out that further exploration is needed to obtain information on the true position of a nearby resonance.
In the two subsequent papers [4] and the more detailed version [7], Rittby et al. combined complex scaling with some Weyl-type analysis and numerical integration methods to compute 44 resonance approximations, including approximations for the first resonance and bound state suggested in [6]; the second resonance therein was further studied by Engdahl & Brändas in [8] by computing lower bounds for norms of Riesz projections. The main conclusion of [4,7] is that there exists a complex threshold thresh with Re( thresh ) > 0, Im( thresh ) < 0 such that all resonance approximations of (1.1) satisfy Re( ) ≤ Re( thresh ) ∼ 4.68 and beyond this threshold, i.e. for Im( ) < Im( thresh ) < 0, their real parts exhibit a certain oscillatory behaviour.
Shortly after the publication of [4] and the submission of [7], Korsch et al. announced in [5, comment] that they had computed a different set of resonance approximations beyond the threshold which did not exhibit any oscillatory behaviour, whereas their earlier computations of lower resonances in [9] had not shown such a disagreement. They used a complex-rotated Milne method and they believed to have backed up their computations by some WKB approximations. Korsch et al. concluded that the results of Rittby et al. for higher resonances were incorrect; they conjectured this might be due to numerical instabilities or to the too limited range 0 < θ < π/4 of angles in the complex scaling method in [4,7].
In an immediate reply (see [10, reply to comment]), Rittby et al. [5] defended their results and attributed the discrepancies of the results to wrongly chosen outgoing boundary conditions. They argued that the asymptotic solutions of the complex Riccati equation associated with (1.1) undergo a dramatic change when θ passes the critical value θ crit = π/4 of the potential in (1.1)  [4,7] against variations of the rotation angle θ, Rittby et al. [4,7] believed to have found approximations to true resonances. About 10 years later, Andersson corroborated the arguments and conclusions of Rittby et al. by a careful multiple-transition point WKB analysis and explained the failure of the complex-rotated Milne method of Korsch et al. by semi-classical theory in [11].
Almost 20 years after the 1982 dispute, the resonance problem (1.1) was studied as an example in two papers in the mathematical literature. In [12], for more general classes of exponentially decaying potentials, Brown et al. developed a resonance-finding procedure for resonances close to points of spectral concentration on the real axis. This method relies on analytic continuation of the Weyl-Titchmarsh function rather than on complex scaling and was first used by Hehenberger et al. [13] in numerical computations for the Stark effect. As an example, Brown et al. computed approximations to the first three resonances of (1.1) which were very close to the ones found in [7]; note that μ/h 2 = 1 in [7] and that the potential q and spectral parameter λ in [12] are related to the potential V and spectral parameter in (1.1) by Not long after, Abramov et al. [14] proved some global analytical bounds for resonances for various classes of potentials. They combined complex scaling with operator theoretic techniques such as numerical ranges and Birman-Schwinger type arguments. Moreover, for the particular case of (1.1), they also performed numerical computations. The analytical results in [14] supported the conjecture of Rittby et al. that a wrong asymptotic boundary condition was used by Korsch et al. [5]. The numerical results of [14] reproduced the resonances found in [4,7] and they suggested three pairs of additional resonances. Each pair consists of an even and an odd resonance so close to each other that they could not be computed accurately. These new resonance pairs may be related to the oscillatory behaviour of the real parts; because two of these pairs satisfy −9.57 ∼ Im( thresh ) < Im( ) < 0.
As it was rightly put in [14], none of the above methods for finding resonances can be used to locate them accurately, but there is clear evidence that they exist. Moreover, none of these methods allows for a proof that a numerically computed candidate for a resonance is not near any true resonance.
The method presented here permits us to settle both questions definitely and adds new information on the threshold beyond which oscillatory behaviour of the real parts of resonances occurs. We prove that the 44 numerical approximations of resonances from [4,7] do lie near true resonances and that the numerical approximations labelled 16-28 in [5] do not lie near true resonances. Moreover, we prove that two of the additional pairs of resonances conjectured in [14] do exist. Our provably correct computations are based on a combination of two key tools, the argument principle on the analytic side and interval arithmetic on the computational side.
Briefly, our approach is as follows. By means of complex scaling x → e iθ x with θ ∈ [0, π/4), the resonances of (1.1) are given in terms of the eigenvalues z = e 2iθ (2 − 1.6) of a Sturm-Liouville problem on R with complex potential. These eigenvalues can be characterized as the zeros of an analytic function . Hence, their number in a rectangle R 0 can be counted by means of the argument principle. On the other hand, we can compute the contour integral in the argument principle in interval arithmetic, using a code based on the software library VNODE developed by Nedialkov et al. [15]. Roughly speaking, this means that all computations, from adding numbers up to integration, amount to working with two-sided estimates; e.g. the sum of two real numbers a ∈ [a 1 , a 2 ] and b ∈ [b 1 , b 2 ] is the interval [a 1 + b 1 , a 2 + b 2 ] which is guaranteed to contain a + b (see [16, §2] for a more detailed description). If we obtain that

Complex scaling and lack of analytic information
There are various mathematical definitions of resonances and different methods to study them; for details, we refer to the comprehensive review articles by Simon [17], Siedentop [18] and Harrell [19]. Here, we use the method of complex scaling where resonances are characterized as eigenvalues of certain non-self-adjoint Schrödinger operators.
As an example, we consider the spectral problem (1.1), with μ/h 2 = 1 for the sake of simplicity. If we set λ := 2 − 1.6, it is easy to see that (1.1) is equivalent to the spectral problem for the linear operator L in the Hilbert space L 2 (R) given by note that W 2 2 (R) is the closure of C ∞ 0 (R) with respect to the norm of W 2 2 (R) given by y 2,2 := ( y 2 2 + y 2 2 + y 2 2 ) 1/2 , where y , y denote the weak derivatives and · 2 denotes the norm of L 2 (R) ( [20, ch. V]).
Because the potential q θ is complex-valued and hence all the above operators H θ along with H D θ , H N θ are no longer self-adjoint, numerical approximations of eigenvalues-and hence of resonances-are prone to be unstable. Examples for such instabilities may be found in [23] for resonances, but they occur already for eigenvalues of matrices (see e.g. [16] for the famous Godunov matrix).
Thus, it is sufficient to show that every resonance 0 ∈ C with Re( 0 ) > 0.8, Im( 0 ) ≤ 0 belongs to C. For every θ ∈ [0, π/4], the point λ 0 := 2 0 − 1.6 lies in the sector −π/2 < arg λ ≤ 0 and is an eigenvalue of the operatorH θ := e −2iθ H θ given by Because the numerical range of a linear operator contains all eigenvalues, we obtain If we note that q θ (R) = q θ ([0, ∞)) and, in addition to a + (θ), we define then it is easy to see that, for θ ∈ [0, π/4), In particular, every resonance λ 0 of L with −π/2 < arg λ ≤ 0 satisfies {λ ∈ C : Re(λ) sin(2θ) + Im(λ) cos(2θ) ≤ a + (θ)}. Figure 1 shows that the only available analytic information is much too coarse to judge the validity or non-validity of resonance approximations. Therefore, it is necessary to employ a method yielding both guaranteed and much more accurate enclosures and exclosures for eigenvalues of non-self-adjoint eigenvalue problems.

Eigenvalue enclosures for complex-valued potentials
The algorithm we use to establish guaranteed eigenvalue enclosures was developed and described in detail in [16,25]. Briefly, it consists of the following two steps. For the sake of simplicity, we consider the Dirichlet problem (2.4); the approach to the Neumann problem (2.5) is completely analogous.
Step A. Solving a truncated problem with guaranteed error bounds. In order to truncate problem (2.4), we restrict the potential q θ to an interval [0, X] and set it equal to 0 on (X, ∞). The unique (up to scalar multiplication) solution of −y = zy in L 2 ([X, ∞)) is exp(− √ −zx) for Re √ −z > 0. Hence, the problem on [0, X], we have to solve is and y(0) = 0, y (X) = − √ −zy(X). (3.1) The eigenvalues of this regular boundary value problem can be characterized as the zeros of an analytic function and may thus be counted and found by means of the argument principle. The algorithms for the calculation of the analytic function and for the contour integral over a chosen starting box R 0 ⊂ C are performed in interval arithmetic, i.e. with guaranteed error bounds. Having achieved (1.2), we obtain a box that contains a certain number n 0 of eigenvalues of the truncated problem (3.1). Repeating this procedure by suitably subdividing the box R 0 , we may finally arrive at a box R Z of desired precision ε Z that contains exactly one eigenvalue z trunc .
then y 2 (0, z) ∈ [y 2 (0, z)]. By means of the interval arithmetic argument principle already used in step A, we now obtain enclosures for the zeros of y 2 (0, z), and hence for the eigenvalues z true of (2.4) of desired precision. For the above-described method, several parameters have to be provided; in particular, the length X of the truncated interval has to be determined such that α X,θ < 1. To this end, we note that and that [26, 7.1.13] Integrating by parts and substituting t = √ ax, we obtain, for a ≥ 0, Hence, for all X ∈ (0, ∞), θ ∈ [0, π/4) with cos(2θ ) ≥ 0.8/x 2 , we can estimate x 2 e − cos(2θ)x 2 /10 dx ≤ 10 cos(2θ ) e − cos(2θ)X 2 /10 X + 5 cos(2θ)X =: A X,θ (3.3) and we use the analytic expression A X,θ to obtain a rigorous computable upper bound A 0 X,θ for A X,θ and hence for α X,θ , To this end, we first expand cos(2θ ) and use Taylor's theorem with remainder in Lagrange form to see that, for every m ∈ N, If θ is a decimal fraction whose fractional part has three digits, the sum T X,θ (m) is rational and can be evaluated exactly. We choose a rigorous computable lower bound T 0 X,θ (m) of T X,θ (m) as the unique decimal number whose fractional part has six digits and T 0 X,θ (m) + 10 −6 > T X,θ (m) ≥ T 0 X,θ (m) (table 1). The function f (t) := (10/t) e −tX 2 /10 (X + 5/tX), t ∈ (0, 1), is decreasing, hence, again by Taylor's theorem with remainder in Lagrange form, we obtain that, for all m, n ∈ N, Now, we fix m, n ∈ N and proceed in the same way as for T X,θ (m) to obtain a rigorous computable  A 0 X,θ 0 (m, n) for θ 0 := 0.75 < π/4 is also an upper bound of A X,θ (m, n) for θ ∈ (0, θ 0 ). Only in two of our computations (for the resonances numbered 37 − and 44 + ), we needed parameter values θ that are larger than θ 0 = 0.75; their upper bound A 0 X,θ (m, n) is computed separately. We use X = 50, m = 2, n = 32 and obtain the rigorous computable lower bounds T 0 X,θ (m) =: T 0 X,θ and upper bounds A 0 X,θ (m, n) =: A 0 X,θ displayed in table 1; note that for X = 50 the condition cos(2θ) ≥ 0.8/x 2 allows for θ ≤ 0.5 arccos( 8 25 10 −5 ), e.g. θ ≤ 0.785238 very close to π/4 ∼ 0.7853981635.

Guaranteed resonance enclosures and exclosures
In [10, reply to comment], Rittby et al. listed a set of 44 approximate resonances ± k of (1.1) that they computed numerically, along with a set of 40 approximate resonances claimed to be found numerically by Korsch et al. in [5, comment]; here, the superscript + occurs for even k, whereas − occurs for odd k. The differences in modulus between these two approximate resonance strings are smaller than 2 · 10 −3 up to − 15 and start to be larger than 10 −2 from + 16 on, getting as huge as 56. 19 [5, comment]. In addition, we enclosed the two pairs of resonances discovered numerically in [14] that are visible by the complex scaling method.
All computed enclosures for resonances, except for one of these pairs, were performed with interval length X = 50, varying scaling angle θ as displayed in the tables, and corresponding guaranteed upper bound A 0 X,θ for α X,θ as in table 1 at the end of §3. The enclosure for one of the additional resonance pairs in [14] turned out to be by far more challenging than all other computations.
We computed guaranteed enclosures for the two pairs of resonances nearλ 2 ,λ 3 ; each of these pairs originates in one eigenvalue z of (2.4) with boundary condition y(0) = 0 with odd eigenfunction (denoted by superscript '−') and one eigenvalue z of (2.5) with boundary condition     for each of the two boundary conditions there is only one resonanceλ − 3 andλ + 3 , respectively, in the disjoint boxes displayed in table 6. Moreover, they guarantee that in the larger λ-box e −2iθ ([23, 24] + [0.05, 1]i) containing these two boxes as well as the numerical valueλ 3 of Abramov et al. there is only one resonance for each of the two boundary conditions. Altogether, we thus proved that there is precisely one pair of disjoint resonancesλ − 3 =λ + 3 near the resonance approximationλ 3 = 2.46 − 23.22i of Abramov et al. and that this approximation has distance approximately 1 · 10 −1 to the true resonance pairλ ± 3 . The computation of the resonance pairλ − 2 ,λ + 2 turned out to be much harder and computationally more expensive than all other enclosures and exclosures. To make it work, we had to use a slight modification of usual complex scaling, using stretching by some parameter R > 0 in addition to rotation of the variable by an angle θ ∈ [0, π/4). The potential q θ,R and the eigenvalue parameter z in the spectral problem for the corresponding operator H θ,R (compare (2.2), (2.3)) then become q θ,R (x) := R 2 e 2iθ (R 2 e 2iθ x 2 − 1.6) e −R 2 e 2iθ x 2 /10 , x ∈ R, z := R 2 e 2iθ λ; note that usual complex scaling corresponds to R = 1.
The main benefit of the additional stretching is that the upper bound A X,θ,R decays exponentially fast in R. As for usual complex scaling, we then applied Taylor's theorem with remainder in Lagrange form to obtain the rigorous computable upper bound A 0 X,θ,R = 1.77 · 10 −17 for X = 10, θ = 0.76 and R = 10.
Hence, our guaranteed enclosures prove that not far from each of the two numerically computed valuesλ 2 andλ 3 of Abramov et al. there is indeed a pair of true resonances of (2.1); the distance is approximately 2 · 10 −2 forλ 2 and approximately 1 · 10 −1 forλ 3 .

Conclusion
In this paper, we have presented a method which, for the first time, permits one to compute resonances in atomic physics with absolute certainty. At the same time, it allows one to detect with absolute certainty wrongly computed resonance approximations. The absolute reliability of our approach is based on a combination of interval arithmetic and the argument principle. To prove the efficiency of our method, we have established guaranteed enclosures for all numerical resonance approximations of Rittby et al. in [4,7] for problem (1.1) and guaranteed exclosures for the numerically computed values of Korsch et al. in [5] that are visible to complex scaling, thus definitely settling a dispute between these two groups of authors. The greatest challenge was to provably enclose two additional pairs of approximate resonances computed by Abramov et al. in [14]