text:
\begin{document}

\maketitle

\begin{abstract}
The question of whether it is possible to compute scattering resonances of Schr\"odinger operators -- independently of the particular  potential -- is addressed. A positive answer is given, and it is shown that the only information required to be known \emph{a priori} is the size of the support of the potential. The potential itself is merely required to be $\mathcal{C}^1$. The proof is constructive, providing a universal algorithm which only needs to access the values of the potential at any requested point.
\end{abstract}

\setcounter{tocdepth}{1}

\section{Introduction and Main Result}
This paper provides an affirmative answer to the following question:
	\begin{quote}\emph{Does there exist a  universal algorithm for computing the  resonances of Schr\"odinger operators with complex potentials?}
	\end{quote}

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 \emph{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. \cite{Hansen11,AHS,Ben-Artzi2015a}) and draws inspiration from   the seminal result \cite{DM} on solving quintic equations via a \emph{tower of algorithms}. {\bf We therefore emphasize that ours is an abstract result in analysis,  not  in numerical analysis.}

\subsection{Quantum Resonances}\label{sec:background}
Let us first define what a quantum resonance is. Let $q:\R^d\to\C$ be compactly supported, let
	\begin{equation*}
	H_q:=-\Delta+q
	\end{equation*}
be the associated Schr\"odinger operator in $L^2(\R^d)$ and let $\chi:\R^d\to\R$ be some compactly supported function with $\chi\equiv 1$ on $\supp(q)$.
It follows from the explicit form of the free fundamental solution (cf. eq. \eqref{eq:Fundamental_Solution} below) that the map
	\begin{equation*}
	z\mapsto I + q(-\Delta-z^2)^{-1}\chi
	\end{equation*}
is an analytic operator-valued function on $\C\setminus\{0\}$. We define:

\begin{de}[Resonance]\label{def:res}
 A \emph{resonance} of $H_q$ is defined
to be a pole of the meromorphic operator-valued function $z\mapsto (I + q(-\Delta-z^2)^{-1}\chi)^{-1}$.
\end{de}
This definition is independent of the specific choice of $\chi$ (so long as $\chi\equiv 1$ on $\supp(q)$), and coincides with the poles of the scattering matrix of $q$, cf. \cite[Prop. 8]{Melrose} and \cite[III.5]{LP}.\\

Resonances can be regarded as states whose wave function disperses very slowly in time, and can therefore be considered as ``almost bound states''.
In physics, such phenomena arise in the description of unstable particles and radioactive decay. Resonant states, just like eigenfunctions, can only exist at  certain energies. The slow-dispersal-in-time approach to resonances motivates one of the earlier definitions of resonances used in the computational physics literature, namely maximization of the so-called time delay function -- see, e.g., Le Roy and Liu \cite{LRL} and Smith \cite{Smith}. This approach leads to real resonance energies for real-valued potentials and, in the one-dimensional case at least, is closely related to the concept of spectral concentration -- see, e.g., Eastham \cite{MSPE}, which describes one mechanism by which such concentrations may arise. For additional discussion we refer to the review article \cite{Zworski1999} and the book \cite{DZ}.

It is widely accepted that the  reliable computation of resonances is a challenging task. This is not usually due to the intrinsic ill-posedness of analytic continuation, since that step is usually done explicitly, but rather due to the fact that complex scaling changes resonance problems either into non-selfadjoint spectral problems, for which the pseudospectra may be far from the spectrum \cite{ET},  or into problems with a nonlinear dependence on the spectral parameter, for which sensitivity to perturbations may also be problematic.  In this context we refer to \cite{Bogli2014}  (including the references and discussion therein) where interval-arithmetic was used to compute resonances.

We show that resonances can be computed as the limit of a sequence of approximations, each of which can be computed precisely using finitely many arithmetic operations. The proof is constructive: we define an algorithm and prove its convergence. We emphasize that this \emph{single} algorithm is valid for \emph{any} Schr\"odinger operator $H_q$ as defined above, so long as   $\diam(\supp(q))$  satisfies an \emph{a priori} bound. We  implement this algorithm in one-dimension and compare its performance to that of Bindel and Zworski \cite{Bindel}.

 \subsection{The Solvability Complexity Index}\label{subsec:sci}
The Solvability Complexity Index (SCI) addresses questions which are at the nexus of pure and applied mathematics, as well as computer science:
 	\begin{quote}\emph{How do we compute objects that are ``infinite'' in nature if we can only handle a finite amount of information and perform finitely many mathematical operations? Indeed, what do we even mean by ``computing'' such an object?}
	\end{quote}
These broad topics are addressed in the sequence of papers \cite{Hansen11,AHS,Ben-Artzi2015a}. Let us summarize the main definitions and discuss how these relate to our problem of finding resonances:
\begin{de}[Computational problem]\label{def:computational_problem}
	A \emph{computational problem} is a quadruple $(\Om,\Lambda,\Xi,\mathcal M)$, where
	\begin{enumi}
		\item $\Om$ is a set, called the \emph{primary set},
		\item $\Lambda$ is a set of complex-valued functions on $\Om$, called the \emph{evaluation set},
		\item $\mathcal M$ is a metric space,
		\item $\Xi:\Om\to \mathcal M$ is a map, called the \emph{problem function}.
	\end{enumi}
\end{de}

\begin{de}[Arithmetic algorithm]\label{def:Algorithm}
	Let $(\Om,\Lambda,\Xi,\mathcal M)$ be a computational problem. An \emph{arithmetic algorithm} is a map $\Gamma:\Om\to\mathcal M$ such that for each $T\in\Om$ there exists a finite subset $\Lambda_\Gamma(T)\subset\Lambda$ such that
	\begin{enumi}
		\item the action of $\Gamma$ on $T$ depends only on $\{f(T)\}_{f\in\Lambda_\Gamma(T)}$,
		\item for every $S\in\Om$ with $f(T)=f(S)$ for all $f\in\Lambda_\Gamma(T)$ one has $\Lambda_\Gamma(S)=\Lambda_\Gamma(T)$,
		\item the action of $\Gamma$ on $T$ consists of performing only finitely many arithmetic operations on $\{f(T)\}_{f\in\Lambda_\Gamma(T)}$.
	\end{enumi}
\end{de}

\begin{de}[Tower of arithmetic algorithms]\label{def:Tower}
	Let $(\Om,\Lambda,\Xi,\mathcal M)$ be a computational problem. A \emph{tower of algorithms} of height $k$ for $\Xi$ is a family $\Gamma_{n_1,n_2,\dots,n_k}:\Om\to\mathcal M$ of arithmetic algorithms such that for all $T\in\Om$
	\begin{align*}
		\Xi(T) = \lim_{n_k\to\infty}\cdots\lim_{n_1\to\infty}\Gamma_{n_1,n_2,\dots,n_k}(T).
	\end{align*}
\end{de}

\begin{de}[SCI]
	A computational problem $(\Om,\Lambda,\Xi,\mathcal M)$ is said to have a \emph{Solvability Complexity Index ($\SCI$)} of $k$ if $k$ is the smallest integer for which there exists a tower of algorithms of height $k$ for $\Xi$.
	If a computational problem has solvability complexity index $k$, we write \begin{align*}
 			\SCI(\Om,\Lambda,\Xi,\mathcal M)=k.
		 \end{align*}
\end{de}

In the present article, our computational problem is made up of the following elements (to be specified more precisely in Section \ref{sec:Setting}):
	\begin{enumi}
		\item $\Om$ is a class of Schr\"odinger operators $H_q$ with potentials $q$ which have a common compact support and a uniform bound in $\mathcal C^1$,
		\item $\Lambda$ is the set of all pointwise evaluations of $q$, as well as pointwise evaluations of the Green's function associated with the Helmholtz operator $-\Delta-z^2$,
		\item $\mathcal M$ is the space $\mathrm{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),
		\item $\Xi:\Om\to \mathcal M$ is the map that associates to a particular Schr\"odinger operator its set of resonances, and we denote it by $\Res(H_q)$.
	\end{enumi}

We show that for this computational problem there exists a tower of height $1$, i.e. there exists a family of algorithms $\{\Gamma_n\}_{n\in\N}$ such that $\Gamma_n(H_q)\to\Res(H_q)$ as $n\to+\infty$ for any $H_q\in\Omega$, where the convergence is in the sense of the Attouch-Wets metric \cite{beer}, generated by the following distance function:

\begin{de}[Attouch-Wets distance] Let $A,B$ be closed sets in $\C$. The  \emph{Attouch-Wets distance} between them is defined as
\begin{align*}
	d_{\text{AW}}(A,B) = \sum_{i=1}^\infty 2^{-i}\min\left\{ 1\,,\,\sup_{|x|<i}\left| \dist(x,A) - \dist(x,B) \right| \right\}.
\end{align*}
Note that if $A,B\subset\C$ are bounded, then $d_{\text{AW}}$ is equivalent to the Hausdorff distance.
\end{de}
\begin{remark}\label{remark:Attouch-Wets}
	It can be shown  (cf. \cite[Prop. 2.8]{R19}) that a sequence of sets $A_n\subset \C$ converges to $A$ in Attouch-Wets metric, if the following two conditions are satisfied
\begin{itemize}
	\item If $\lambda_n\in A_n$ and $\lambda_n\to\lambda$, then $\lambda\in A$.
	\item If $\lambda\in A$, then there exist $\lambda_n\in A_n$ with $\lambda_n\to\lambda$.
\end{itemize}

\end{remark}

\subsection{Main Result}\label{sec:Setting}

Let $d\in\mathbb N$, fix $M>0$ and let $\Om_M$ denote the class of Schr\"odinger operators
\begin{align*}
	H_q:=-\Delta+q\quad\text{on}\quad L^2(\R^d)
\end{align*}
with $q\in \mathcal C^1(\R^d;\C)$ and $\supp(q)\Subset Q_M$, where $Q_M$ denotes the cube of edge length $M$ centered at the origin.
Moreover, for $x\in\R^d$ and $z\in\C$ let
\begin{align}\label{eq:Fundamental_Solution}
	G(x,z):=\begin{cases}\f\i 4\left( \f{z}{2\pi |x|} \right)^{\f{d-2}{2}}H^{(1)}_{\f{d-2}{2}}\big( z|x| \big), &d\ge 2,\\
		\f{\i}{2z}e^{\i z|x|}, & d=1,
		\end{cases}
\end{align}
where $H^{(1)}_{\nu}$ denotes the Hankel function of the first kind. For $\im(z)>0$, $G(x,z)$ is the fundamental solution to the free Helmholtz operator $-\Delta-z^2$ (cf. \cite[Ch. 22]{Serov}).
We define the evaluation set $\Lambda$ to be
\begin{align}\label{eq:input_info}
	\Lambda &:=  \left\{M \right\}\cup\left\{ q\mapsto q(x)\,|\, x\in\R^d \right\}
		\cup \left\{  G(x,z)\,|\, x\in\R^d,\,z\in\C \right\}.
\end{align}
Then the quadruple $(\Om_M,\Lambda,\Res(\cdot),\text{cl}(\C))$ poses a computational problem in the sense of Definition \ref{def:computational_problem}. The main result of the present article is the following.
\begin{theorem}\label{th:mainth}
Resonances can be computed in one limit: $\SCI(\Om_M,\Lambda,\Res(\cdot),\mathrm{cl}(\C))=1.$
\end{theorem}
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 \ref{sec:Numerics}.

\begin{remark}
\begin{enumi}
	\item We note that there exist examples of computational spectral problems for which $\SCI\geq2$, even in the selfadjoint case  (cf. \cite[Th. 6.5]{AHS}). Whenever $\SCI\geq2$ 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.
	\item If we drop the restriction $\supp(q)\Subset Q_M$ from the definition of $\Omega$, i.e. if the size of $\supp(q)$ is not assumed to be known \emph{a priori}, then Theorem \ref{th:mainth} implies $\SCI\leq 2$, but it is no longer clear that  $\SCI =1$. Indeed, let us provide a sketch of the proof:

		\emph{Sketch of proof:}  Choose a sequence $0\leq M_k\to+\infty$ and for each fixed $k$ let $\{\Gamma_{k,n}\}_{n\in\N}$ be a sequence of algorithms as in Theorem \ref{th:mainth}, which computes the resonances of $\Om_{M_k}$ in one limit. Given a compactly supported potential function $q$, whose support is not known, choose a smooth cutoff function $\rho_k$ with $Q_{M_k-1}\subset\supp(\rho_k)\Subset Q_{M_k}$. Then from Theorem \ref{th:mainth} we know that $\{\Gamma_{k,n}\}_{n\in\N}$ computes the set of resonances of $H_{\rho_k q}$ in one limit. As soon as $\supp(q)\Subset Q_{M_k-1}$, one has $\rho_k q\equiv q$ and the sequence $\{\lim_{n\to\infty}\Gamma_{nk}(H_{\rho_kq})\}_{k\in\N}$ will be constant in $k$, and hence convergent. This proves that $\lim_{k\to\infty}\lim_{n\to\infty}\Gamma_{k,n}(q)=\Res(H_q)$, that is, the set of resonances can be computed in two limits.
\end{enumi}
\end{remark}

The proof of Theorem \ref{th:mainth} is divided into several steps. First, we obtain quantitative resolvent norm estimates for the operator $K(z):=q(H_q-z^2)^{-1}\chi$ from Definition \ref{def:res}. 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 $(I+K_n(z))^{-1}$.

\subsubsection*{Organization of the paper}

Section \ref{sec:Analytic_Continuation} contains a short discussion of Definition \ref{def:res} and meromorphic continuation. In Section \ref{sec:abstract} we prove some estimates for convergence of finite-dimensional approximations of linear operators, which are then used in Section \ref{sec:algorithm} to construct an explicit algorithm which computes resonances in one limit, thereby proving Theorem \ref{th:mainth}. Section  \ref{sec:Numerics} is dedicated to numerical experiments. In Appendix \ref{app:greens} we review some properties of the Green's function $G(x,z)$ introduced in \eqref{eq:Fundamental_Solution}.

\section{Analytic Continuation}\label{sec:Analytic_Continuation}
We  use this section for a more detailed discussion of Definition \ref{def:res} and to fix some notations and conventions.
For the sake of self containedness, we  prove the existence of $z\mapsto (I+q(-\Delta-z^2)^{-1}\chi)^{-1}$ as a meromorphic operator-valued function on the domain
\begin{align*}
	\C^{\mathrm{ext}}:=\begin{cases}
		\C & \text{if } d \text{ is odd},\\
		\text{logarithmic cover of }\C & \text{if } d \text{ is even}.
	\end{cases}
\end{align*}
This result follows from the classical Analytic Fredholm Theorem (cf. e.g. \cite[Sec. VI.5]{RS})
\begin{theorem}[Analytic Fredholm Theorem]
	Let $D\subset \C$ be open and connected and let $F:D\to L(\h)$ be an analytic operator-valued function such that $F(z)$ is compact for all $z\in D$. Then, either
	\begin{enumi}
		\item $(I+F(z))^{-1}$ exists for no $z\in D$, or
		\item $(I+F(z))^{-1}$ exists for all $z\in D\setminus S$, where $S$ is a discrete subset of $D$. In this case, $z\mapsto (I+F(z))^{-1}$ is meromorphic in $D$, analytic in $D\setminus S$, the residues at the poles are finite rank operators, and if $z\in S$ then $\ker(I+F(z))\neq\{0\}$.
	\end{enumi}
\end{theorem}

Next, recall that $Q_M$ denotes the cube of edge length $M$ in $\R^d$ centered at the origin. Let $\chi:=\chi_{Q_M}$ be the indicator function of $Q_M$. Note that the operator-valued function $z\mapsto q(-\Delta-z^2)^{-1}\chi$ is an analytic function on $\C^{\mathrm{ext}}\setminus\{0\}$. This follows from the explicit representation of the free fundamental solution \eqref{eq:Fundamental_Solution} (cf. Remark \ref{remark:even_odd}).

\begin{lemma}
	The function $\C^+\ni z\mapsto\big(I+q(-\Delta-z^2)^{-1}\chi\big)^{-1}$ has a meromorphic continuation to $\C^{\mathrm{ext}}$.
	Moreover, the residues at the poles are finite rank operators.
\end{lemma}
\begin{proof}
The operator $q(-\Delta-z^2)^{-1}\chi$ is compact by the Fr\'echet-Kolmogorov theorem and the inverse $\big(I+q(-\Delta-z^2)^{-1}\chi\big)^{-1}$ exists for $\im(z)>0$ large enough, by the Neumann series. Hence, the claim follows from the analytic Fredholm theorem, together with Remark \ref{remark:even_odd} in the appendix.
\end{proof}

The above observations lead us to study the spectrum of the compact operator
\begin{align}\label{eq:free_resolvent}
	K(z):=q(-\Delta-z^2)^{-1}\chi, \qquad z\in\mathbb \C^{\mathrm{ext}}.
\end{align}
Since the integral kernel for the free resolvent is given explicitly by \eqref{eq:Fundamental_Solution} as an analytic function of $z\in\C^{\text{ext}}\setminus\{0\}$, we have an explicit representation of \eqref{eq:free_resolvent} as an integral operator on $L^2(\R^d)$:
\begin{align}\label{eq:integral_formula}
	\big(q(-\Delta-z^2)^{-1}\chi\, f\big)(x) = q(x)\int_{\R^d} G(x-y,z)\chi(y)f(y)\,dy, \qquad z\in \C^{\mathrm{ext}}\setminus\{0\}.
\end{align}

\section{Abstract Error Estimates}\label{sec:abstract}
We recall that  the resonances of $H_q=-\Delta+q$ are defined to be  the poles of $\C^{\mathrm{ext}}\ni z\mapsto\big(I+K(z)\big)^{-1}$ where $K(z)=q(-\Delta-z^2)^{-1}\chi$ is a compact operator. 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 \ref{th:mainth}. 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\subset \h$ be a finite-dimensional subspace, $P_n:\h\to\h_n$ the orthogonal projection  and $K:\C^{\mathrm{ext}}\to L(\h)$ continuous in operator norm.  Moreover, let $K_n:\C^{\mathrm{ext}}\to L(\h_n)$ be analytic for every $n\in\N$. Assume that for any compact subset $B\subset \C^\mathrm{ext}$ there exist a sequence $a_n\downarrow0$ and a constant $C>0$ such that
\begin{align}
	\|K(z)-K_n(z)P_n\|_{L(\h)} &\leq Ca_n, \label{eq:K-KnPn}\\
	\left\|P_nK(z)|_{\h_n}-K_n(z)\right\|_{L(\h_n)} &\leq Ca_n, \label{eq:PnK|-Kn}\\
	\|K(z)-P_nK(z)P_n\|_{L(\h)} &\leq Ca_n,\label{eq:K-PKP}
\end{align}
for all $z\in B$.

\subsection{Error Estimates}
\begin{lemma}\label{lemma:pollution}
	If $z\in\C^\mathrm{ext}$ is such that $-1\notin\sigma(K(z))$, then

	\begin{align*}
		\big(1-Ca_n\|(I+K(z))^{-1}\|_{L(\h)}\big)\left\|(I+K_n(z))^{-1}\right\|_{L(\h_n)} \leq \left\|(I+K(z))^{-1}\right\|_{L(\h)},
	\end{align*}
	where we use the convention that $\|(I+K_n(z))^{-1}\|_{L(\h)}=+\infty$ if $-1\in\sigma(K_n(z))$.
\end{lemma}
\begin{proof}
	Whenever the left hand side is non-positive the assertion is trivially true, so we may assume w.l.o.g. that $1-Ca_n\|(I+K(z))^{-1}\|_{L(\h)}>0$.
	In this case, the assertion follows by a Neumann series argument, as follows. We have
	\begin{equation}\label{eq:factor}
	\begin{aligned}
		I+K_n(z)P_n &= I+K(z)+(K_n(z)P_n-K(z))\\
		&= (I+K(z))\left[ I+(I+K(z))^{-1}(K_n(z)P_n-K(z)) \right]
	\end{aligned}
	\end{equation}
	Because $Ca_n<\f{1}{\|(I+K(z))^{-1}\|}$, the second factor in \eqref{eq:factor} is invertible by the Neumann series and
	\begin{align*}
		\left[ I+(I+K(z))^{-1}(K_n(z)P_n-K(z)) \right]^{-1} = \sum_{j=0}^\infty \Big((I+K(z))^{-1}(K_n(z)P_n-K(z))\Big)^j.
	\end{align*}
	Hence,
	\begin{align*}
		\left\|(I+K_n(z)P_n)^{-1}\right\|_{L(\h)} &\leq \Bigg\| \sum_{j=0}^\infty \Big((I+K(z))^{-1}(K_n(z)-K(z))\Big)^j \Bigg\|_{L(\h)} \left\|(I+K(z))^{-1}\right\|_{L(\h)} \\
		&\leq \sum_{j=0}^\infty \left\|(I+K(z))^{-1}\right\|_{L(\h)}^{j+1} \left\|K_n(z)P_n-K(z)\right\|_{L(\h)}^j\\
		&\leq \sum_{j=0}^\infty \left\|(I+K(z))^{-1}\right\|_{L(\h)}^{j+1} (Ca_n)^j\\
		&= \left\|(I+K(z))^{-1}\right\|_{L(\h)} \sum_{j=0}^\infty \left\|(I+K(z))^{-1}\right\|_{L(\h)}^{j} (Ca_n)^j\\
		&= \f{\|(I+K(z))^{-1}\|_{L(\h)}}{1-\left\|(I+K(z))^{-1}\right\|_{L(\h)} Ca_n}
	\end{align*}
for any $n\in\N$.	It remains to replace the $L(\h)$ norm on the left hand side by the $L(\h_n)$ norm. This follows from Claim \ref{claim:replace-norm}. 	This completes the proof.
\end{proof}
\begin{cl}\label{claim:replace-norm}
	 We have $\|(I+K_n(z))^{-1}\|_{L(\h_n)}\leq \|(I+K_n(z)P_n)^{-1}\|_{L(\h)}$ for all $z$ for which both operators are boundedly invertible.
\end{cl}
	\begin{proof}
	For $x\in\h_n$ we have $(I+K_nP_n)^{-1}x = (I+K_n)^{-1}x$, because if $u\in\h_n$ solves $(I+K_n)u=x$, then $(I+K_nP_n)u=x$ and by invertibility it follows that $u=(I+K_nP_n)^{-1}x$. We conclude that
		\begin{align*}
			\sup_{x\in\h_n,\|x\|=1} \|(I+K_nP_n)^{-1}x\|_{\h} = \sup_{x\in\h_n,\|x\|=1} \|(I+K_n)^{-1}x\|_{\h_n}
		\end{align*}
		and therefore
		\begin{align*}
			\sup_{x\in\h,\|x\|=1} \|(I+K_nP_n)^{-1}x\|_{\h} \ge \sup_{x\in\h_n,\|x\|=1} \|(I+K_n)^{-1}x\|_{\h_n}.
		\end{align*}
\end{proof}

\begin{lemma}\label{lemma:inclusion}
	If $z\in\mathbb C^\mathrm{ext}$ is such that either $-1\in\sigma(K(z))$ or $\left\|(I+K(z))^{-1}\right\|_{L(\h)}\ge \f{1}{Ca_n}$, then either $-1\in\sigma(P_nK(z)P_n)$ or
	\begin{equation*}\label{eq:resolvent>1/3a}
		\left\|(I+P_nK(z)P_n)^{-1}\right\|_{L(\h)}\ge \f{1}{2Ca_n}.
	\end{equation*}
\end{lemma}
\begin{proof}
	If $-1\in\sigma(K(z))$, then unless $-1\in\sigma(P_nK(z)P_n)$, we have
	\begin{equation*}\label{eq:factor2}
	\begin{aligned}
		I+K(z) &= I+P_nK(z)P_n+(K(z)-P_nK(z)P_n)\\
		&= (I+P_nK(z)P_n)\left[ I+(I+P_nK(z)P_n)^{-1}(K(z)-P_nK(z)P_n) \right]
	\end{aligned}
	\end{equation*}
	We now argue by contradiction. If we had $\|(I+P_nK(z)P_n)^{-1}\|_{L(\h)}<\f{1}{2Ca_n}$, then we would have $\|(I+P_nK(z)P_n)^{-1}(K(z)-P_nK(z)P_n)\|_{L(\h)}<1$ and $I+K(z)$ would be invertible by the Neumann series contradicting our assumption that $-1\in\sigma(K(z))$. Thus we must have $\|(I+P_nK(z)P_n)^{-1}\|_{L(\h)} \geq \f{1}{2Ca_n}$.

	Now let us turn to the case where $-1\notin\sigma(K(z))$ and $\left\|(I+K(z))^{-1}\right\|_{L(\h)}\ge \f{1}{Ca_n}$. The same calculation as in the proof of Lemma \ref{lemma:pollution} shows that
	\begin{align*}
		\left(1-Ca_n\left\|(I+P_nK(z)P_n)^{-1}\right\|_{L(\h)}\right)\left\| (I+K(z))^{-1} \right\|_{L(\h)} \leq \left\|(I+P_nK(z)P_n)^{-1}\right\|_{L(\h)}
	\end{align*}
	from which it follows easily that $\f{1}{2Ca_n}\leq \left\|(I+P_nK(z)P_n)^{-1}\right\|_{L(\h)}$.
\end{proof}

\subsection{An Abstract Algorithm For Computing Poles}

We now demonstrate how the  the assumptions \eqref{eq:K-KnPn}-\eqref{eq:K-PKP} allow us to construct an abstract algorithm that computes the poles of  $\big(I+K(z)\big)^{-1}$. By an \emph{abstract algorithm} we mean a sequence of subsets of $\C^{\mathrm{ext}}$, which is constructed from $K_n$ and which converges in Attouch-Wets metric to $\{z\in \C^{\mathrm{ext}}\,|\,-1\in\sigma(K(z))\}$. Note that this is not yet an \emph{arithmetic algorithm} in the sense of Definition \ref{def:Algorithm}, since the sets are not computed from a finite amount of information in finitely many steps.

Let $B\subset\C^{\mathrm{ext}}$ be compact and define the exponentially fine lattice $\mathcal{L}_n:=e^{-\f1{a_n}}(\Z+\i\Z)\cap 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
\begin{align*}
	\Theta_n^B(K) = \left\{z\in \mathcal{L}_n\,\bigg|\,\left\|(I+K_n(z))^{-1}\right\|_{L(\h_n)} \ge \f{1}{2\sqrt{a_n}} \right\}.
\end{align*}
Moreover, note that by \cite[Prop. 10.1]{AHS}, determining whether $\left\|(I+K_n(z))^{-1}\right\|_{L(\h_n)} \ge \f{1}{2\sqrt{a_n}}$ can be done with finitely many arithmetic operations on the matrix elements of $K_n(z)$ for each $z\in\mathcal L_n$.
\begin{lemma}\label{lemma:convergence}
	The assumptions \eqref{eq:K-KnPn}-\eqref{eq:K-PKP} imply the convergence $\Theta_n^B(K)\to \{z\in B\,|\,-1\in\sigma(K(z))\}$ in Attouch-Wets metric.
\end{lemma}

\begin{proof}

\emph{I. Excluding spectral pollution.} Assume that $z_n\in \Theta_n^B(K)$ with $z_n\to z_0$ for some $z_0\in B$. Then for each $n$ we have $\|(I+K_n(z_n))^{-1}\|_{L(\h_n)} \ge \f{1}{2\sqrt{a_n}}$ and hence by Lemma \ref{lemma:pollution}
\begin{equation*}
\left\|(I+K(z_n))^{-1}\right\|_{L(\h)} \ge \Big(1-Ca_n\left\|(I+K(z_n))^{-1}\right\|_{L(\h)}\Big)\f12a_n^{-\f12}.
\end{equation*}
(with the convention that $\|(I+K(z_n))^{-1}\|_{L(\h)}=+\infty$ if $-1\in\sigma(K(z_n))$). Whenever $\sqrt{a_n}\leq \f{2}{C}$ this leads to
\begin{align*}
	\left\|(I+K(z_n))^{-1}\right\|_{L(\h)} &\ge \f12\, \f{a_n^{-\f12}}{1+C\f{\sqrt{a_n}}{2}}\\
	&\ge \f14 a_n^{-\f12}.
\end{align*}
It follows  that $\|(I+K(z_n))^{-1}\|_{L(\h)} \to+\infty$ as $n\to+\infty$ 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.

\emph{II. Spectral inclusion.}
Assume now that $z$ is a pole, i.e. $-1\in\sigma(K(z))$. Our reasoning will have the structure
\begin{align*}
	-1\in & \sigma(K(z))\\
	&\Downarrow\\
	\exists z_n\in \mathcal{L}_n : \|(I&+K(z_n))^{-1}\|_{L(\h)} \;\text{ large} \\
	&\Downarrow\\
	\|(I+P_nK&(z_n) P_n)^{-1}\|_{L(\h)} \;\text{ large}\\
	&\Downarrow\\
	\|(I+P_nK&(z_n)|_{\h_n})^{-1}\|_{L(\h_n)} \;\text{ large}\\
	&\Downarrow\\
	\|(I+K_n&(z_n))^{-1}\|_{L(\h_n)} \;\text{ large,}
\end{align*}
with a quantitative estimate in each step. To this end, note first that if $-1\in\sigma(K(z))$ for some $z\in B$, then there exist $\nu,c,\eps>0$ (independent of $n$) such that for all  $\zeta$ in a $\eps$-neighborhood of $z$,
\begin{align}\label{eq:Laurent}
	\|(I+K(\zeta))^{-1}\|_{L(\h)}\geq c |z-\zeta|^{-\nu}.
\end{align}
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 \eqref{eq:Laurent} that for any $z_n$ such that $|z-z_n|\leq e^{-\f1{a_n}}$ one will have
\begin{equation*}
	\|(I+K(z_n))^{-1}\|_{L(\h)} \geq c |z-z_n|^{-\nu}
	\geq c e^{\f{\nu}{a_n}}
	\geq \f1{Ca_n}
\end{equation*}
for all large enough $n$. We conclude that for any pole $z$ there exists a sequence $z_n\in \mathcal{L}_n$ such that  $z_n\to z$ as $n\to+\infty$ and $\|(I+K(z_n))^{-1}\|_{L(\h)}>\f1{Ca_n}$ for  all but finitely many $n\in\N$.

Next, Lemma \ref{lemma:inclusion} shows that $\|(I+P_nK(z_n)P_n)^{-1}\|_{L(\h)}>\f{1}{2Ca_n}$. Studying this norm further, we have
\begin{align*}
	(I_\h+P_nK(z_n)P_n)^{-1} &= \bigl(I_{\h_n}+P_nK(z_n)|_{\h_n}\bigr)^{-1}\oplus I_{\h_n^\perp}
\end{align*}
and thus
\begin{align*}
	\left\|(I_\h+P_nK(z_n)P_n)^{-1}\right\|_{L(\h)} &= \max\left\{\bigl\|\bigl(I_{\h_n}+P_nK(z_n)|_{\h_n}\bigr)^{-1}\bigr\|_{L(\h_n)}\,,\, 1\right\}.
\end{align*}
Hence, as soon as $a_n<\f1{2C}$, we have $\|(I+P_nK(z_n)P_n)^{-1}\|_{L(\h)} = \|(I+P_nK(z_n)|_{\h_n})^{-1}\|_{L(\h_n)}$. We conclude that if $z$ is a pole, then there exists $z_n\in \mathcal{L}_n$ such that
\begin{align}\label{eq:PnK|_bound}
	\|(I+P_nK(z_n)|_{\h_n})^{-1}\|_{L(\h_n)}>\f1{2Ca_n}
\end{align}
($n$ large enough). A similar reasoning as in Lemma \ref{lemma:pollution} (using \eqref{eq:PnK|-Kn}) shows that now
\begin{align*}
	\big(1-Ca_n\|(I+K_n(z_n))^{-1}\|_{L(\h_n)}\big)\|(I+P_nK(z_n)|_{\h_n})^{-1}\|_{L(\h_n)} \leq \|(I+K_n(z_n))^{-1}\|_{L(\h_n)},
\end{align*}
and rearranging terms, together with \eqref{eq:PnK|_bound}, gives
\begin{align*}
	\|(I+K_n(z_n))^{-1}\|_{L(\h_n)} \geq \f{1}{4Ca_n}
\end{align*}
and therefore $z_n\in\Theta_n^B(K)$ for large enough $n$.
The assertion about Attouch-Wets convergence now follows from {Remark \ref{remark:Attouch-Wets}.}
\end{proof}

\section{Definition of the Algorithm}\label{sec:algorithm}

\subsection{Error Estimates}

In this section, we will apply the abstract results of Section \ref{sec:abstract} to our resonance problem. To this end, define $K(z) = q(-\Delta-z^2)^{-1}\chi$.
We write $K$ for the integral kernel of $K(z)$ to simplify notation. Recall that $K$ was given by \eqref{eq:integral_formula} and $\supp(K)\subset Q_M\times Q_M$.
We will construct an operator approximation $K_n$ of $K$, which satisfies \eqref{eq:K-KnPn}-\eqref{eq:K-PKP} and in addition
\begin{itemize}
	\item[(H1)] The matrix elements of $K_n$ can be computed in finitely many steps from a finite subset $\Lambda_n\subset\Lambda$ (cf. eq. \eqref{eq:input_info} and Def. \ref{def:Algorithm});
	\item[(H2)] The convergence rate $a_n$ is explicitly known (i.e. the sequence $a_n$ can be used to define the algorithm).
\end{itemize}
To this end, let us define $\h_n,P_n$ as follows.
\begin{align}
	\R^d &= \bigcup_{i\in\f1n\Z^d} S_{n,i} :=\bigcup_{i\in\f1n\Z^d} \left(\big[0,\tfrac 1 n\big)^d+i\right) \label{eq:R^d_decomp}\\
	\h_n &= \left\{ f\in L^2(Q_M)\,\middle|\, f|_{S_{n,i}} \text{ constant }\forall i\in\tfrac1n\Z^d \cap Q_M \right\}\notag \\
	P_nf(x) &= \sum_{i\in\f1n\Z^d\cap Q_M} \bigg(n^d\int_{S_{n,i}} f(t)\,dt\bigg)\chi_{S_{n,i}}(x)\label{eq:P_n_def}
\end{align}
Furthermore, we have to make a concrete choice for the approximation $K_n$. An obvious choice is the integral kernel
\begin{align*}
	K_n(x,y) = \sum_{i,j\in n^{-1}\Z^d\cap Q_M} K(i,j)\chi_{S_{n,i}}(x)\chi_{S_{n,j}}(y),
\end{align*}
i.e. a piecewise constant approximation of $K(\cdot,\cdot)$ 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. \eqref{eq:K-KnPn}-\eqref{eq:K-PKP}. To streamline the presentation, we will restrict ourselves to $d\ge 3$ in our computations, the cases $d\leq 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.

\paragraph{Proof of \eqref{eq:K-PKP}.}

Using the definitions \eqref{eq:R^d_decomp}-\eqref{eq:P_n_def}, we have
\begin{align*}
	Kf(x) - P_nKP_nf(x) &= \int_{\R^d} K(x,y)f(y)\,dy - \int_{\R^d} P_n^xK(x,y)P_nf(y)\,dy,
\end{align*}
where $P_n^xK(x,y)$ means $(P_nK(\cdot,y))(x)$. Using $L^2$-selfadjointness of $P_n$, we conclude that
\begin{align*}
	Kf(x) - P_nKP_nf(x) &= \int_{\R^d} K(x,y)f(y)\,dy - \int_{\R^d} P_n^yP_n^xK(x,y)f(y)\,dy\\
	 &= \int_{\R^d} \big(K(x,y)- P_n^yP_n^xK(x,y)\big)f(y)\,dy,
\end{align*}

Note that $P_n^yP_n^xK(x,y)$ simply yields a step function approximation of $K(x,y)$ like \eqref{eq:P_n_def}, but in dimension $2d$. We conclude by applying Young's inequality \cite[Th. 0.3.1]{Sogge}, that
\begin{align*}
	\|Kf - P_nKP_nf\|_{L^2(\R^d)} &\leq \eta_n \|f\|_{L^2(\R^d)},
\end{align*}
where
\begin{align}\label{eq:eta_n}
	\eta_n = \max\left\{ \sup_{x\in\R^d}\int_{\R^d} |K(x,y)- P_n^yP_n^xK(x,y)|\,dy \,,\,\sup_{y\in\R^d}\int_{\R^d} |K(x,y)- P_n^yP_n^xK(x,y)|\,dx \right\}
\end{align}
Thus, all we have to do is estimate the $L^\infty$-$L^1$ difference between $K$ and its projection onto step functions. To this end, fix $x\in Q_M$, let $\eps>\f{2}{n}$ and decompose the integrals as follows
\begin{align}\label{eq:K-PPK_whole}
\begin{aligned}
	\int_{\R^d} |K(x,y)- & P_n^yP_n^xK(x,y)|\,dy = \int_{Q_M} |K(x,y)-P_n^yP_n^xK(x,y)|\,dy  \\
	&= \int_{Q_M\setminus B_\eps(x)} |K(x,y)-P_n^yP_n^xK(x,y)|\,dy + \int_{B_\eps(x)} |K(x,y)-P_n^yP_n^xK(x,y)|\,dy
\end{aligned}
\end{align}
The integral over $B_\eps(x)$ can be estimated by $\int_{B_\eps(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\in\f1n\Z^d$ be such that $x\in S_{n,j}$, see Figure \ref{fig:Sketch}. Let  $i\in\f1n\Z^d$ be such that $|i-j|>\f{\eps}{2}$. Then:

\begin{figure}[htbp]
	\centering
	\input{Integral_Sketch.tex}
	\caption{
		Sketch of the geometry in the calculation leading to \eqref{eq:complement_of_B_eps}. The sum over $i$ includes all cells whose nodes are outside the dashed ball centered at $j$.
	}
	\label{fig:Sketch}
\end{figure}

\begin{align*}
	\int_{S_{n,i}} |K(x,y) & -P_n^yP_n^xK(x,y)|\,dy = \int_{S_{n,i}} \left|K(x,y)-\dashint_{S_{n,i}\times S_{n,j}} K(s,t)dsdt\right|\,dy\\
	&\leq \int_{S_{n,i}} \dashint_{S_{n,i}\times S_{n,j}} \left|K(x,y)- K(s,t)\right|dsdt\,dy\\
	&= \int_{S_{n,i}} \dashint_{S_{n,i}\times S_{n,j}} \left|\int_0^1\nabla K\big(\tau\big(\begin{smallmatrix}x\\y\end{smallmatrix}\big)+(1-\tau)\big(\begin{smallmatrix}s\\t\end{smallmatrix}\big)\big)\cdot\big(\big(\begin{smallmatrix}x\\y\end{smallmatrix}\big)-\big(\begin{smallmatrix}s\\t\end{smallmatrix}\big)\big)d\tau \right|dsdt\,dy\\
	&\leq \int_{S_{n,i}} \dashint_{S_{n,i}\times S_{n,j}} \int_0^1\big|\nabla K\big(\tau\big(\begin{smallmatrix}x\\y\end{smallmatrix}\big)+(1-\tau)\big(\begin{smallmatrix}s\\t\end{smallmatrix}\big)\big)\big|\big|\big(\begin{smallmatrix}x\\y\end{smallmatrix}\big)-\big(\begin{smallmatrix}s\\t\end{smallmatrix}\big)\big| d\tau\, dsdt\,dy\\
	&\leq \int_{S_{n,i}} \dashint_{S_{n,i}\times S_{n,j}} \int_0^1\big\|\nabla K\big\|_{L^\infty(S_{n,i}\times S_{n,j})}\f{2\sqrt d}n d\tau\, dsdt\,dy\\
\end{align*}
Summing over $i$, we finally obtain (cf. Figure \ref{fig:Sketch})
\begin{align}
	\int_{\R^d\setminus B_{\eps}(x)} |K(x,y)-P_n^yP_n^x K(x,y)|\,dy &\leq \sum_{i:|i-j|>\f\eps2} \int_{S_{n,i}} |K(x,y)-P_n^yP_n^xK(x,y)|\,dy\nonumber\\
	&\leq \sum_{i:|i-j|>\f\eps2} \int_{S_{n,i}} \dashint_{S_{n,i}\times S_{n,j}} \int_0^1\big\|\nabla K\big\|_{L^\infty(S_{n,i}\times S_{n,j})}\f{2\sqrt d}n d\tau\, dsdt\,dy\nonumber\\
	&\leq \f{2\sqrt d}n \big\|\nabla K\big\|_{L^\infty(Q_M\setminus B_{\f{\eps}{2}}(x))} \int_{Q_M\setminus B_{\f{\eps}{4}}(x)} \,dy\nonumber\\
	&= |Q_M|\f{2\sqrt d}{n} \big\|\nabla K\big\|_{L^\infty(Q_M\setminus B_{\f{\eps}{2}}(x))} \nonumber\\
	&\leq  |Q_M|\f{2\sqrt d}{n} \|q\|_{\mathcal C^1} C\left(\f{\eps}{2}\right)^{1-d}\nonumber\\
	&\leq C\f{|Q_M|}n \eps^{1-d}, \label{eq:complement_of_B_eps}
\end{align}
where the fifth line follows from \eqref{eq:G_bound_2}, in the appendix, and the bound $\|q\|_{\mathcal C^1}\leq +\infty$. Using \eqref{eq:complement_of_B_eps} in \eqref{eq:K-PPK_whole}, we conclude that
\begin{align*}
	\int_{\R^d} |K(x,y)- P_n^yP_n^xK(x,y)|\,dy &\leq C\f{|Q_M|}{n} \eps^{1-d} + \int_{B_\eps(x)} 2|K(x,y)|\,dy\\
	\Rightarrow\quad\sup_{x\in \R^d}\int_{\R^d} |K(x,y)- P_n^yP_n^xK(x,y)|\,dy &\leq C\f{|Q_M|}{n} \eps^{1-d} + C'\eps^2,
\end{align*}
where in the last line we have used \eqref{eq:G_bound} and the boundedness of $q$ again.

With an analogous calculation for $\sup_{y\in\R^d}\int_{\R^d} |K(x,y)- P_n^yP_n^xK(x,y)|\,dx$ (which we omit here), and recalling that $\eta_n$ was defined by \eqref{eq:eta_n}, we conclude that for all $\eps>0$
\begin{align*}
	\eta_n\leq \f{1}{n}C\eps^{1-d} + C'\eps^2.
\end{align*}
Choosing $\eps:=n^{-\f{1}{d+1}}$, we conclude that
\begin{align}\label{eq:explicit_rate}
	\|Kf - P_nKP_nf\|_{L^2(\R^d)} &\leq \f{C+C'}{n^{\f{2}{d+1}}} \|f\|_{L^2(\R^d)}
\end{align}
and hence $\|K-P_nKP_n\|_{L(L^2(\R^d))}\to 0$ as $n\to+\infty$ with rate (at least) $a_n=n^{-\f{2}{d+1}}\leq n^{-\f1d}$.

\begin{remark}
	Note that the constants $C,C'$ all depend on the spectral parameter $z$, but are bounded for $z$ in compact subsets of $\C^{\mathrm{ext}}$, because $K$ depends continuously on $z$.
\end{remark}

\paragraph{Proof of \eqref{eq:PnK|-Kn} and (H1).}

An orthonormal basis of $\h_n$ is given by the functions
\begin{align*}
	e_i:=n^{\f d2}\chi_{S_{n,i}},\qquad i\in\tfrac1n\Z^d\cap Q_M,
\end{align*}
so that
	\begin{equation*}
	P_nf = \sum_{j\in \f1n\Z\cap Q_M}\p{f,e_j}_{L^2}e_j
	\end{equation*}
in this basis. It is then easily seen that in this basis $K_n$ has the matrix elements
\begin{equation*}
	(K_n)_{ij} = n^{-d} K(i,j).
\end{equation*}
Note that this proves (H1): The matrix elements of $K_n$ can be calculated in finitely many arithmetic operations from the finite set $\Lambda_n:=\{K(i,j)\,|\,i,j\in\f1n\Z\cap Q_M\}\subset\Lambda$.
Similarly, it can be seen that the matrix elements of $P_nK|_{\h_n}$ in this basis are given by
\begin{align*}
	(P_nK)_{ij} &= n^{d}\int_{S_{n,i}}\int_{S_{n,j}}K(x,y)\,dxdy  \\[1mm]
	&=: n^{-d} \langle K\rangle_{ij},
\end{align*}
where we have introduced the notation $\langle \cdot\rangle_{ij}$ for the mean value on $S_{n,i}\times S_{n,j}$.
Let $f=\sum_{j} f_je_j \in \h_n$. From the above, and Young's inequality, we conclude that
\begin{align*}
	\|(P_nK-K_n)f\|_{L^2}^2 &= \sum_{i\in \f1n\Z^d\cap Q_M}\left|\sum_{j\in \f1n\Z^d\cap Q_M} n^{-d} \big(K(i,j)-\langle K\rangle_{ij}\big)f_j\right|^2\\
	&\leq \tilde\eta_n^2\|f\|_{L^2}^2,
\end{align*}
where
\begin{align*}
		\tilde\eta_n := \max\Bigg\{ \sup_{i\in\f1n\Z^d\cap Q_M}\,\sum_{j\in \f1n\Z^d\cap Q_M} n^{-d}|K(i,j)-\langle K\rangle_{ij}| \,,\,\sup_{j\in\f1n\Z^d\cap Q_M}\,\sum_{i\in \f1n\Z^d\cap Q_M} n^{-d}|K(i,j)-\langle K\rangle_{ij}| \Bigg\}.
\end{align*}
Hence, we have reduced the problem to estimating these $\ell^\infty$-$\ell^1$ differences. This can be done similarly to \eqref{eq:K-PPK_whole}, by separating $(Q_M\times Q_M)\cap (\f1n\Z\times\f1n\Z)$ into an $\eps$-region around $i=j$ and the rest:
\begin{align}
	\sum_{j\in \f1n\Z^d\cap Q_M} n^{-d}|K(i,j)-\langle K\rangle_{ij}| &= \sum_{|j-i|>\eps}n^{-d}|K(i,j)-\langle K\rangle_{ij}| + \sum_{|j-i|\leq\eps}n^{-d}|K(i,j)-\langle K\rangle_{ij}| \nonumber \\
	&\leq Cn^{-1}\sum_{|j-i|>\eps} n^{-d} \|\nabla K\|_{L^\infty(\{|x-y|>\eps\})} + \sum_{|j-i|\leq\eps}n^{-d}|K(i,j)-\langle K\rangle_{ij}|  \nonumber\\
	&\leq Cn^{-1} \eps^{-d+1} + \sum_{|j-i|\leq\eps}n^{-d}|K(i,j)-\langle K\rangle_{ij}|,  \label{eq:discrete_estimate}
\end{align}
where we have used \eqref{eq:G_bound_2} and the $\mathcal C^1$-boundedness of $q$ in the last line. To estimate the last term on the right hand side, note that $|K(i,j)-\langle K\rangle_{ij}|\leq C|j-i|^{-(d-2)}$ near $i=j$ (cf. eq. \eqref{eq:G_bound}).
Next, note that the sum $n^{-d}\sum_{j:|j-i|\leq\eps} \f{1}{|j-i|^{d-2}}$ can be interpreted as an integral over a piecewise constant function, which approximates $(x,y)\mapsto |x-y|^{2-d}$. But this function is dominated by $(x,y)\mapsto |x-y|^{1-d}$ when $|x-y|$ is small, and therefore we have
\begin{align}
	n^{-d}\sum_{j:|j-i|\leq\eps} \f{1}{|j-i|^{d-2}} &\leq C \int_{B_{2\eps}(x)} |x-y|^{1-d}\,dy \nonumber\\
	&= C\int_0^{2\eps}  r^{1-d}\,\omega_d r^{d-1} dr \nonumber\\
	&= 2C\omega_d\, \eps \label{eq:discrete_2}
\end{align}
where $\omega_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 \eqref{eq:discrete_2} into \eqref{eq:discrete_estimate}, we arrive at
\begin{align*}
	\sum_{j\in \f1n\Z^d\cap Q_M} n^{-d}|K(i,j)-\langle K\rangle_{ij}|
	&\leq Cn^{-1} \eps^{-d+1} + 2C\omega_d\, \eps.
\end{align*}
Choosing $\eps=n^{-\f1d}$ yields
\begin{align}\label{eq:explicit_rate2}
	\sum_{j\in \f1n\Z^d\cap Q_M} n^{-d}|K(i,j)-\langle K\rangle_{ij}|
	&\leq C'n^{-\f1d}.
\end{align}
 Finally, swapping $i$ and $j$ will give an analogous estimate and we can conclude that $\tilde\eta_n\to 0$ with rate $a_n = n^{-\f1d}$.
\begin{remark}
	Note again that the constants $C,C'$ depend on $z$, but are bounded for $z$ in compact subsets of $\C^{\mathrm{ext}}$, since $K$ depends continuously on $z$.
\end{remark}

\paragraph{Proof of \eqref{eq:K-KnPn} and (H2).}
Estimate \eqref{eq:K-KnPn} in fact follows from \eqref{eq:K-PKP} and \eqref{eq:PnK|-Kn}. Indeed, writing $K_n$ and $K$ as block operator matrices w.r.t. the decomposition $\h = \h_n\oplus\h_n^\perp$, we have
\begin{align*}
	K = \begin{pmatrix}
		P_nK|_{\h_n} & D_1\\
		D_2 & D_3
	\end{pmatrix},
\end{align*}
with some operators $D_1,D_2,D_3$. Estimate \eqref{eq:K-PKP} shows that
\begin{align}\label{eq:0BCD}
	\left\|\begin{pmatrix}
		0 & D_1\\
		D_2 & D_3
	\end{pmatrix}\right\|_{L(\h)}<Ca_n,
\end{align}
whereas estimate \eqref{eq:PnK|-Kn} shows that
\begin{align}\label{eq:K000}
	\left\|P_nK|_{\h_n}-K_n\right\|_{L(\h_n)} = \left\|\begin{pmatrix}
		P_nK|_{\h_n}-K_n & 0\\
		0 & 0
	\end{pmatrix}\right\|_{L(\h)}<Ca_n.
\end{align}
Together, eqs. \eqref{eq:0BCD} and \eqref{eq:K000} imply that
\begin{align*}
	\|K(z)-K_n(z)P_n\|_{L(\h)}  = \left\|\begin{pmatrix}
		P_nK|_{\h_n}-K_n & D_1\\
		D_2 & D_3
	\end{pmatrix}\right\|_{L(\h)} < 2Ca_n.
\end{align*}
 The explicit rates obtained in \eqref{eq:explicit_rate} and \eqref{eq:explicit_rate2} prove that our approximation scheme satisfies (H2).

\subsection{The Algorithm}
It remains to extend the algorithm $\Theta_n^B$ from a single compact set $B\subset\mathbb C^\mathrm{ext}$ to the entire complex plane. This is done via a diagonal-type argument.
\subsubsection{Odd Dimensions}
We choose a tiling of $\C$, where we start with a square $B_1=\big\{z\in\C\,\big|\,|\re(z)|\leq\f12, \, -1\leq|\im(z)|\leq 0 \big\}$ and then add squares in a counterclockwise spiral manner as shown in Figure \ref{fig:tiling}.

\begin{figure}[htbp]
	\centering
	\contourlength{2pt}
	\begin{tikzpicture}[>=stealth]
		\draw[-stealth] (-3,0)--(3,0) node[above]{$\mathrm{Re}\, z$};
		\draw[-stealth] (0,-2.5)--(0,2) node[right]{$\mathrm{Im}\, z$};
		\draw   (-0.05,-1)--(0.05,-1);
		\draw   (-0.05,-2)--(0.05,-2);
		\draw   (-0.05,-3)--(0.05,-3);
		\foreach \y in {0,-1,-2}
		\foreach \x in {0,1} {
			\draw (-0.5+\x,-2-\y) rectangle (0.5+\x,-1-\y);
		}
		\draw (-1.5,0) rectangle (-0.5,1);
		\node at (0,-0.5) {\contour{white}{$B_1$}};
		\node at (0,-1.5) {\contour{white}{$B_2$}};
		\node at (1,-1.5) {$B_3$};
		\node at (1,-0.5) {$B_4$};
		\node at (1,0.5) {$B_5$};
		\node at (0,0.5) {\contour{white}{$B_6$}};
		\node at (-1,0.5) {$B_7$};
		\contourlength{1pt}
		\node at (-1,-0.4) {\contour{white}{$\vdots$}};
	\end{tikzpicture}
	\caption{Tiling of the complex plane}
	\label{fig:tiling}
\end{figure}
Next, we define our algorithm as follows. We let
\begin{align*}
	\Gamma_1(q) &:= \Theta_1^{B_1}(q) \\
	\Gamma_2(q) &:= \Theta_2^{B_1}(q) \cup \Theta_2^{B_2}(q) \\
	\Gamma_3(q) &:= \Theta_3^{B_1}(q) \cup \Theta_3^{B_2}(q) \cup \Theta_3^{B_3}(q) \\
	&\vdots\\
	\Gamma_n(q) &:= \bigcup_{j=1}^n \Theta_n^{B_j}(q).
\end{align*}
Lemma \ref{lemma:convergence} ensures that each $\Theta_n^{B_k}$ converges to $\Res(q)\cap B_k$ for fixed $k$ and since the $\{B_k\}$ form a tiling of $\C$, it follows that $\Gamma_n(q)\to \Res(q)$ in Attouch-Wets metric.

\subsubsection{Even Dimensions}

\begin{figure}[htbp]
	\centering
	\contourlength{1pt}
	\input{even_tiling}
	\caption{Tiling of the logarithmic Riemann surface}
	\label{fig:even-dim_tiling}
\end{figure}

In even dimensions we have to cover not only the complex plane $\C$, but its logarithmic covering space, which is equivalent to covering infinitely many copies of the complex plane. A similar strategy as in the odd dimensional case, together with a diagonal-type argument does the job in this case. Indeed, we can construct a cover by boxes $B_n$ as follows (cf. Figure \ref{fig:even-dim_tiling}).
\begin{enumerate}
	\item Start with box $B_1$ (defined as in the odd dimensional case) on the first Riemann sheet;
	\item Add a box $B_2$ below $B_1$ on sheet number 1 and add a box $B_1$ on sheet number 2;
	\item Add a box $B_3$ on sheet number 1, add a box $B_2$ on sheet number 2 and a box $B_1$ on sheet number 3;
	\item \dots
\end{enumerate}

Next, define again
\begin{align*}
	\Gamma_1(q) &:= \Theta_1^{B_1^{(1)}}(q) \\
	\Gamma_2(q) &:= \Theta_2^{B_1^{(1)}}(q) \cup \Theta_2^{B_2^{(1)}}(q) \cup \Theta_2^{B_1^{(2)}}(q) \\
	\Gamma_3(q) &:= \Theta_3^{B_1^{(1)}}(q) \cup \Theta_3^{B_2^{(1)}}(q) \cup \Theta_3^{B_3^{(1)}}(q) \cup \Theta_3^{B_1^{(2)}}(q) \cup \Theta_3^{B_2^{(2)}}(q) \cup \Theta_3^{B_3^{(1)}}(q) \\
	&\vdots\\
	\Gamma_n(q) &:= \bigcup_{k=1}^n\bigcup_{j=1}^{n-k+1} \Theta_n^{B_j^{(k)}}(q).
\end{align*}
Lemma \ref{lemma:convergence} ensures that each $\Theta_n^{B_j^{(k)}}$ converges to $\Res(q)\cap B_j^{(k)}$ for fixed $k$ and since the $\{B_j^{(k)}\}$ form a tiling of $\C^{\mathrm{ext}}$, it follows that $\Gamma_n(q)\to \Res(q)$ in Attouch-Wets metric. The proof of Theorem \ref{th:mainth} is complete.

\section{Numerical Results}\label{sec:Numerics}
Software to compute resonances has been in existence for decades \cite{Rittby81, BEM00, AAD01}. The authors of \cite{BZ} recently proposed a collection of MATLAB codes to compute resonance poles and scattering of plane waves efficiently (``MatScat'', cf. \cite{Bindel}). 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)\subset[a,b]$ (for some known $a<b$), which computes the set
$$\Big\{z\in \mathcal{L}_n\cap B\,\Big|\,\Big\|\Big(\mathds{1}_{n\times n}+\big(K(i,j)\big)_{i,j\in \f{b-a}{n}\Z\cap[a,b]}\Big)^{-1}\Big\|>C\Big\},$$
where the region $B$ in the complex plane, the lattice distance of $\mathcal{L}_n$ and cutoff threshold $C$ were treated as independent parameters.

\paragraph{Comparison of results.}
Figures \ref{fig:Gauss_well} and \ref{fig:Trapping} 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.

\begin{figure}[htbp]
\centering
	\input{Smooth_Gauss}
	\caption{Comparison of the result of \cite{Bindel} (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\times 4M = 1000\times 4000$.}
	\label{fig:Gauss_well}
\end{figure}

\begin{figure}[htbp]
\centering
	\input{Smooth_Trapping}
	\caption{Comparison of the result of \cite{Bindel} (black) and our algorithm (blue) for a smooth trapping potential supported between $-1.2$ and $1.2$. 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\times 10M = 1000\times 10000$.}
	\label{fig:Trapping}
\end{figure}

\paragraph{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 \ref{fig:Gauss_well} 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 there is considerable room for improvement in numerical efficiency.

Moreover, our algorithm can only yield reliable results in a certain region, as the following heuristic calculations make clear.
\begin{itemize}
\item\emph{Imaginary part of $z$:} Since the fundamental solution $G(x,z)=\f{1}{2iz}e^{\i z|x|}$ grows exponentially with $-\im(z)$ and $x\in[-a,a]$, a limit is reached when $|\im(z)|\sim \f{\log(2M)}{2a}$, where $M$ is the largest number the machine can store with adequate precision (for the interval $[-a,a]=[-1,1]$ and $M=10^{16}$ this bound yields $\im(z)\gtrsim -18.8$).

\item\emph{Real part of $z$:} Similarly, a natural bound on $\re(z)$ is reached when the period of $e^{\i z|x|}$ is less than twice the lattice spacing $\f{2}{n}$, i.e. when
$|\re(z)|\lesssim \pi n$ (for $n=30$ this bound yields $|\re(z)|\lesssim 94$).
\end{itemize}
Numerical experiments have confirmed the above bounds (see Figure \ref{fig:Bound_on_real_part}). Note that the bound on $\im(z)$ is fixed by the machine precision, while the bound on $|\re(z)|$ can be raised by increasing $n$.

\begin{remark}
We note that our algorithm is not restricted to one-dimension or real-valued potentials. Indeed,  the algorithm $\Gamma_n$ only uses the bound $\supp(q)\subset Q_M$, and higher dimensional implementations of $\Gamma_n$ can be coded similarly to the one-dimensional one.
\end{remark}
\begin{figure}[htbp]
	\centering
	\input{Bound_Real_Part}
	\caption{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)|\sim \pi n$ in each case.}
	\label{fig:Bound_on_real_part}
\end{figure}

\appendix

\section{Fundamental Solution}\label{app:greens}
In this appendix we  gather some well-known results about the fundamental solution for the Helmholtz equation. These facts are used to show that the abstract framework of Section \ref{sec:abstract} holds in the context of our algorithm as defined in Section \ref{sec:algorithm}, namely that eqs. \eqref{eq:K-KnPn}-\eqref{eq:K-PKP} hold.
We adopt the notation of \cite{AS} and write $f(\zeta)\sim \zeta^\nu$ if $f$ and $\zeta^\nu$ are \emph{asymptotically equal}, i.e. $|f(\zeta)-\zeta^\nu|=O(|\zeta|^{\nu+1})$ as $|\zeta|\to 0$.
\begin{remark}
	By the asymptotic expansion of the Hankel functions
	\begin{equation*}
		H^{(1)}_\nu(\zeta)\sim \begin{cases}
		-\f{\Gamma(\nu)}{\pi}\left(\f{\zeta}{2}\right)^{-\nu}, & \nu>0,\\
 			\f{2\i}{\pi}\log(\zeta), & \nu=0,
	 \end{cases}
	\end{equation*}
	where $\Gamma$ denotes the Gamma function and $\log$ denotes the principal branch of the logarithm (cf. \cite[Ch. 9.1.9]{AS}), we find that the fundamental solution \eqref{eq:Fundamental_Solution} satisfies the small $|x|$ asymptotics
	\begin{align*}
		G(x,z) &\sim -\f{\i\Gamma(\f{d-2}{2})}{\pi}\left(\f{z|x|}{2}\right)^{-\f{d-2}{2}}\f\i 4\left( \f{z}{2\pi |x|} \right)^{\f{d-2}{2}}  \\
		&= \f{\Gamma(\f{d-2}{2})}{4\pi^{\f n 2}}\f{1}{|x|^{d-2}} , \qquad  \text{as }|x|\to 0,
	\end{align*}
	for $n\geq 3$, and
	\begin{align*}
		G(x,z) &\sim -\f{1}{2\pi}\log(z|x|), \qquad \text{as }|x|\to 0,
	\end{align*}
	for $n=2$. Hence
	\begin{align}\label{eq:G_bound}
		|G(x,z)| &\leq C_z\cdot \begin{cases}
		\f{1}{|x|^{d-2}}, &n\geq 3,\\
 		 \log(|x|), &n=2,
	 \end{cases}
	\end{align}
	where $C_z>0$ is uniformly bounded for $z$ in a compact subset of $\C$. Similar formulas hold for the derivatives of $G$. Indeed, identities for Hankel functions (cf. \cite[Ch. 9.1.30]{AS}) show that
	\begin{align}\label{eq:G_bound_2}
		|\nabla G(x,z)| &\leq
		\f{C_z}{|x|^{d-1}},\qquad\text{ for }d\geq 2.
	\end{align}
\end{remark}

\begin{remark}\label{remark:even_odd}
	From the representation of $G(x,z)$ in terms of Hankel functions it follows that $G$ can be  continued analytically in $z$ through the branch cut $\R_+$. In fact, it can be shown that $G$ can be continued to
	\begin{itemize}
		\item the Riemann surface of the complex square root, if $d$ is odd,
		\item the Riemann surface of the complex logarithm, if $d$ is even,
	\end{itemize}
	(cf. \cite[Ch. 3.1.4]{DZ}). The estimates \eqref{eq:G_bound} and \eqref{eq:G_bound_2} remain valid in either case.
\end{remark}


TITLE: Gaussian formula for $n$ dimensions
QUESTION [0 upvotes]: I know this is a very simple one. If this is the formula for the two dimensional Gaussian (no covariance matrix considered - I have one mean and variance for each dimension):
$$ A\exp{\left[ -\left(\displaystyle \frac{(\mu_{1} - x_{1})^2}{2\sigma_{1}^2} + \frac{(\mu_{2} - x_{2})^2}{2\sigma_{2}^2}\right) \right]}$$
Would this be the one for the three dimensional and so on?
$$ A\exp{\left[ -\left(\displaystyle \frac{(\mu_{1} - x_{1})^2}{2\sigma_{1}^2} + \frac{(\mu_{2} - x_{2})^2}{2\sigma_{2}^2} + \frac{(\mu_{3} - x_{3})^2}{2\sigma_{3}^2} + \dots \right) \right]}$$
Thanks.

REPLY [0 votes]: This would indeed be correct:
$$ A\exp{\left[ -\left(\displaystyle \frac{(\mu_{1} - x_{1})^2}{2\sigma_{1}^2} + \frac{(\mu_{2} - x_{2})^2}{2\sigma_{2}^2} + \frac{(\mu_{3} - x_{3})^2}{2\sigma_{3}^2} + \dots \right) \right]}$$
Please have a look at the comments if you want to ensure that things integrate to one.


TITLE: Absolute value: First Derivative Heaviside Function + Second Derivative Dirac Delta Function Distribution
QUESTION [1 upvotes]: I have some issues concerning the derivative of an absolute value $|x|$, the Heaviside function $\theta(x)$ and the Dirac Delta Distribution.
Given the definition of the Heaviside function
\begin{equation}
  \theta(x) = \left\{
                \begin{array}{ll}
                  1 \quad \text{for} \quad x \ge 0\\
                  0 \quad \text{for} \quad x<0
                \end{array}
              \right.
\end{equation}
I suppose it is right to express the first derivative of an absolute value as
\begin{equation}
  \frac{d |x|}{dx} = \left\{
                \begin{array}{ll}
                  +1 \quad \text{for} \quad x > 0\\
                  -1 \quad \text{for} \quad x<0
                \end{array}
              \right.
\overset{!}{=} \theta(x) - \theta(-x),
\end{equation}
evaluating the derivatives of $|x|$ from the left and right side. Here I have some issues with the case of $x=0$
\begin{equation}
  \left.\frac{d |x|}{dx}\right|_{x=0}=\theta(0)-\theta(0)=1-1=0,
\end{equation}
which is not defined in my original "first derivative" of $d |x|/ dx$.
So how can the first derivative of an absolute value be correctly expressed in terms of the Heaviside function?
Anyways taking my assumption of the first derivative for granted I want to perform a second derivative with the identity
\begin{equation}
  \frac{d \theta(x)}{dx} = \delta(x)
\end{equation},
which then leads to
\begin{equation}
  \frac{ d^2 |x|}{dx^2} = \frac{d}{dx}(\theta(x)-\theta(-x)) = \delta(x)-\delta(-x)=\delta(x)-\delta(x)=0.
\end{equation}
The last expression cannot be true and should be
\begin{equation}
  \frac{d^2 |x|}{dx^2} = 2 \delta(x)
\end{equation}
following Second derivative of absolute value function proportional to Dirac delta function?
What am I missing to get to the correct expression?

REPLY [0 votes]: It does make more sense than you think, I think.
Although there is no true derivative for the absolute value, Ive seen optimization problems such as this treated.  I do believe it was called the "subgradient".
Basically you take the one-sided derivatives at $0$, and you get both $+1$ and $-1$ simultaneously for $x=0$, depending on the side you approach from.  This forms a range of values, I will express as the interval $(-1,1)$.  An argument can be made that $f' \in (-1,1)$, being all possible slopes concurrently.  This can be argued because any line you draw through the point at $x=0$ and having any slope in this range will fail to cross through the function.  And as you can see, $0\in(-1,1)$; a slope of $0$ is indicative of a peak. Consequently, $x=0$ is a local optimization point. So its not unreasonable to suggest that the derivative of $|x|$ at $x=0$ is zero, or, at least, we can say that a horizontal line is tangent to the function (i.e. touches it once without crossing) and that a peak exists there.


TITLE: Special Relativity Problem
QUESTION [0 upvotes]: I am having trouble with the following problem:

Fry travels in a rocket ship towards Leela, at constant relative speed $v$: Fry
  is delivering a pizza, which in its rest frame stays hot for exactly another
  2 minutes. If Leela measures that Fry is 27 million kilometers away, then
  calculate the minimal value of $v$ for which the pizza is hot when delivered.

My approach was to use speed=distance/time and account for time dilation but I cannot figure out how to eliminate the Lorentz Factor. The answer is meant to be $3c/5$. Any pointers would be appreciated! Thanks.

REPLY [0 votes]: Since an answer has already been posted using time dilation and $\gamma$, here's an alternative method employing the invariant interval.
$$(c\tau)^2 = (c \Delta t)^2 - \Delta x^2$$
As specified, the proper time for the pizza is $\tau = 120s$
The displacement of the pizza is $\Delta x = 27 \cdot 10^6 km = 90$ light-seconds.
Solving for the elapsed coordinate time $\Delta t$ yields
$$\Delta t = \sqrt{ \tau^2 + \frac{\Delta x^2}{c^2}} = \sqrt{120^2 + 90^2} = 150s$$
Thus,
$$v = \frac{\Delta x}{\Delta t} = \frac{90}{150}c= \frac{3}{5}c$$


TITLE: Nearly Consistent, is the following true?
QUESTION [0 upvotes]: Continuing my previous question:

Let A be a set of axioms, we say that A is nearly consistent is there
is a finite sub set of A (A') for which A\A' is consistent.

Is it true to say that:
A is nearly consistent iff each finite sub set of A is  nearly consistent?
Which side is correct and which one isn't?

REPLY [1 votes]: The answer is essentially contained in the answers/comments to your previous question:

Using the fact that a subset of a consistent set is consistent, we can easily prove that any subset of a nearly consistent set is nearly consistent. (This is basically the content of J.G.'s answer.) So the left-to-right direction holds.

In the other direction, in the comments below my answer I mentioned that every finite set of sentences is nearly consistent (do you see why?). So no matter what $A$ is, every finite subset of $A$ is nearly consistent. From this it's easy to get a counterexample to the right-to-left direction:

 Just take $A$ to be any set which is not nearly consistent.

(Note that the fact that every finite set is nearly consistent makes my argument in the first bulletpoint unnecessary, but it's still worth understanding.)


TITLE: Operators and power regularity
QUESTION [3 upvotes]: Let $E$ be an infinite-dimensional complex Hilbert space.

An operator $A\in \mathcal{L}(E)$ is called power regular if for all $x\in E$, $r(x,A):=\lim_{n\rightarrow \infty} \|A^nx\|^{1/n}$ exits.

Assume that for all $x,y\in E$, $r(x,A),r(y,A)$ and $r(x+y,A)$ exists.

I want to show that
  $$r(x+y,A)\leq \max\{r(x,A),r(y,A)\} .$$

REPLY [2 votes]: We have
\begin{align}
\|A^n(x+y)\|^{1/n} &= \|A^nx+A^ny\|^{1/n} \\
&\le \Big(\|A^nx\|+\|A^ny\|\Big)^{1/n} \\
&\le \Big(\max\{\|A^nx\|, \|A^ny\}\| + \max\{\|A^nx\|, \|A^ny\}\|\Big)^{1/n}\\
&= \Big(2\max\{\|A^nx\|, \|A^ny\}\|\Big)^{1/n}\\
&= 2^{1/n }\max\left\{\|A^nx\|^{1/n}, \|A^ny\|^{1/n}\right\}
\end{align}
Letting $n\to\infty$ gives
\begin{align}
r(x+y, A) &= \lim_{n\to\infty}\|A^n(x+y)\|^{1/n} \\
&\le \lim_{n\to\infty} 2^{1/n }\max\left\{\|A^nx\|^{1/n}, \|A^ny\|^{1/n}\right\} \\
&= \max\left\{\lim_{n\to\infty}\|A^nx\|^{1/n}, \lim_{n\to\infty}\|A^ny\|^{1/n}\right\}\\
&= \max\{r(x, A), r(y, A)\}
\end{align}