Computing the Sound of the Sea in a Seashell

The question of whether there exists an approximation procedure to compute the resonances of any Helmholtz resonator, regardless of its particular shape, is addressed. A positive answer is given, and it is shown that all that one has to assume is that the resonator chamber is bounded and that its boundary is C2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal {C}}}^2$$\end{document}. The proof is constructive, providing a universal algorithm which only needs to access the values of the characteristic function of the chamber at any requested point.


Introduction
This paper provides an affirmative answer to the following question: Does there exist a universal algorithm for computing the resonances of the Laplacian in R 2 \ U for any open bounded set U ⊂ R 2 ?Any domain U ⊂ R d gives rise to resonances, i.e. special frequencies that are 'nearly' eigenvalues of the Laplacian in R d \ U (with appropriate boundary conditions).The most famous example is the "sound of the sea" in a seashell: when we hold a seashell against our ear we hear frequencies that are nearly the eigenvalues of the Laplacian in the closed cavity with Neumann boundary conditions.This is also known as a Helmholtz resonator [18].We therefore refer to U as a "resonator".
We are interested in computing (Dirichlet) resonances in a way that is independent of the domain U itself; i.e.U is the input of the problem, and the computation returns the associated resonances.We work in the plane, i.e. d = 2, but our results can be adapted to higher dimensions.To our 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.[15,4]) 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.Resonances.Let us first define what we mean by a resonance.Let U ⊂ R d be an open set and assume that ∂U ∈ C 2 .Let H be the Laplacian in L 2 (R d \ U ) with homogeneous Dirichlet boundary conditions on ∂U .Resonances of H can be defined via analytic continuation of the associated Dirichlet-to-Neumann (DtN) operators.Indeed, we can start from a characterization of eigenvalues and perform analytic continuation as follows.Denote the spectral parameter by k 2 (with the branch cut of the complex square root running along the positive real line).Then k 2 is an eigenvalue of H if and only if there exists a function u ∈ L 2 (R d \ U ), such that (−∆ − k 2 )u = 0 in R d \ U and u = 0 on ∂U .The existence of such a u is equivalent to the following.Let R > 0 be such that U ⊂ B R (the open ball of radius R around 0). Denote by M inner (k) and M outer (k), respectively, the inner and outer DtN maps on ∂B R associated with −∆ − k 2 .M inner (k) and M outer (k) can be shown to be analytic operator-valued functions of k in the upper half plane C + := {z ∈ C | Im(z) > 0}.If there exists an eigenfunction u of H with eigenvalue k 2 , then u| ∂B R ∈ ker(M inner (k) + M outer (k)).
Conversely, if φ ∈ ker(M inner (k) + M outer (k)), then there exist solutions u inner and u outer of (−∆ − k 2 )u = 0 in B R \U and R d \B R respectively, such that u inner | ∂B R = u outer | ∂B R = φ and ∂ ν u inner | ∂B R = −∂ ν u outer | ∂B R (note that the normal vector ν is always taken to point to the exterior of the domain considered, that is, away from 0 for B R and towards 0 for R d \ B R ).Hence, the function is in H 2 (R d \ U ).This shows that a complex number k ∈ C + in the upper half plane is a (square root of an) eigenvalue if and only if ker(M inner (k) + M outer (k)) = ∅.A resonance can now be defined as follows.Moreover, an argument similar to the one surrounding eq.(1.1) shows that the nontriviality of ker(M inner (k) + M outer (k)) is independent of the value of R as long as U ⊂ B R .
Remark 1.3.In two dimensions, the DtN maps M inner , M outer can actually be analytically continued to the logarithmic cover of C, rather than just C − (cf.[12,Sec. 3.1.4]).However, since the vast majority of the literature is concerned with resonances near the real axis, we decided to simplify our presentation by considering only resonances in C − .A strategy for dealing with the general case has been outlined in [6,Sec. 4.2.2].
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 and accessing finitely many values of the characteristic function χ U of U .The proof is constructive: we define an algorithm and prove its convergence.We emphasize that this single algorithm is valid for any open U ⊂ R 2 with ∂U ∈ C 2 .We implement this algorithm and compare its output to known results.
1.2.The Solvability Complexity Index.The Solvability Complexity Index (SCI) addresses questions which are at the nexus of pure and applied mathematics, as well as computer science: 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?These broad topics are addressed in the sequence of papers [15,4,5].Let us summarize the main definitions and discuss how these relate to our problem of finding resonances: Definition 1.4 (Computational problem).A computational problem is a quadruple (Ω, Λ, Ξ, M), where (i) Ω is a set, called the primary set, (ii) Λ is a set of complex-valued functions on Ω, called the evaluation set, (iii) M is a metric space, (iv) Ξ : Ω → M is a map, called the problem function.
Definition 1.5 (Arithmetic algorithm).Let (Ω, Λ, Ξ, M) be a computational problem.An arithmetic algorithm is a map Γ : Ω → M such that for each T ∈ Ω there exists a finite subset Λ Γ (T ) ⊂ Λ such that (i) the action of Γ on T depends only on {f (T )} f ∈ΛΓ(T ) , (ii) for every S ∈ Ω with f (T ) = f (S) for all f ∈ Λ Γ (T ) one has Λ Γ (S) = Λ Γ (T ), (iii) the action of Γ on T consists of performing only finitely many arithmetic operations on {f (T )} f ∈ΛΓ(T ) .
Definition 1.6 (Tower of arithmetic algorithms).Let (Ω, Λ, Ξ, M) be a computational problem.A tower of algorithms of height k for Ξ is a family Γ n1,n2,...,n k : Ω → M of arithmetic algorithms such that for all T ∈ Ω Ξ(T ) = lim Definition 1.7 (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.
1.3.Setting of the problem and main result.Let us describe the elements comprising our computational problem, followed by our main theorem.
(i) The primary set Ω. We define (ii) The evaluation set Λ.The evaluation set, which describes the data at our algorithm's disposal, is comprised of (all points in) the set U , as well as some basic special functions which we assume we can access at no additional computational cost and any (nonnegative) real number (such as e, π or √ 2): where n denote the Bessel and Hankel functions of the first kind, respectively.Including the characteristic functions U → χ U (x) in Λ means that for every x ∈ R 2 we can test whether x is included in U or not.
(iii) The metric space M. M is the space cl(C) of all closed subsets of C equipped with the Attouch-Wets metric, generated by the following distance function: Definition 1.8 (Attouch-Wets distance).Let A, B be closed sets in C. The Attouch-Wets distance between them is defined as Note that if A, B ⊂ C are bounded, then d AW is equivalent to the Hausdorff distance d H . Furthermore, it can be shown (cf.[3,Ch. 3]) that (iv) The problem function Ξ.The problem function Ξ(U ) = Res(U ) is the map that associates to each U ∈ Ω the set of resonances (as defined in Definition 1.1) of H.
With these at hand, we can state our main theorem: Theorem 1.9.There exists a family of arithmetic algorithms {Γ n } n∈N such that Γ n (U ) → Res(U ) as n → +∞ for any U ∈ Ω, where the convergence is in the sense of the Attouch-Wets metric.That is, SCI(Ω, Λ, Res(•), (cl(C), d AW )) = 1.
Remark 1.10.In Subsection 5.2 we outline a strategy for adapting this result to Neumann boundary conditions, and implement a numerical model consisting of four circular holes: this is relevant for understanding resonances caused by oil-drilling platforms [13], for instance.Theorem 1.9 is proved by identifying an operator of the form Id L 2 + K(k), where K(k) is in some Schatten class C p , with the property that ker(Id The former is equivalent to det p (Id L 2 + K(k)) = 0, where det p denotes the p -perturbation determinant.This determinant is approximated via a finite element procedure on B R \ U with explicit error bounds.A thresholding procedure then yields an approximation for the zero set of det p (Id L 2 + K(k)).
1.4.Discussion.The study of cavity resonances is much older than the study of quantum mechanical resonances.The foundational work is generally ascribed to Helmholtz [18], who in the 1850s had constructed devices which were designed to identify special frequencies from within a sound wave.These Helmholtz resonators consisted of a small aperture at one end to admit sound, and a larger opening at the other to emit it.
The operator theoretic foundations of the study of cavity resonances are usually based on various approaches to the theory of scattering by obstacles.A good early treatment of both quantummechanical and obstacle scattering may be found in the seminal text of Newton [25].The book of Lax and Phillips [20] is perhaps the most famous work on scattering theory in a semi-abstract setting, though the whole subject was extensively researched in the Soviet school by mathematicians such as Faddeev and Pavlov, see, e.g.[14], and also the monograph of Yafaev [31].
There is extensive work in the applied mathematics literature too, devoted to estimating resonances in various asymptotic regimes, either geometrical or semiclassical: see [16,17,23,30,28] as starting points.Assemblies of resonators are proposed as cloaking devices in a variety of contexts, see, e.g., [1]; the mathematical treatment of acoustic waveguides is very similar to that of optical or quantum waveguides.A very up-to-date overview of the current state of the art in the subject may be found in [12].
Numerical approaches use a variety of techniques of which complex scaling (rediscovered as 'perfectly matched layers'), boundary integral techniques and combinations of these with special finite element methods are the most common.
Organization of the paper .In Section 2 we analyze the inner and outer DtN maps and obtain a new expression equivalent to (1.2) in which M inner (k) + M outer (k) is replaced by an expression of the form Id L 2 +perturbation.Section 3 is dedicated to the construction of a finite element method for the approximation of this perturbation, and in Section 4 the algorithm for approximating Res(U ) is defined and shown to converge.Finally, in Section 5 we provide some numerical examples.

Formulas for the inner and outer DtN maps
As a first step, we prove a weaker version of Theorem 1.9.We show that under the additional assumption that the diameter of the resonator is known a priori, there exists an algorithm that computes the set of resonances in one limit.To this end we define where B ρ is the open ball about the origin of radius ρ > 0. Then one has: Theorem 2.1.There exists a family of arithmetic algorithms {Γ R n } n∈N depending on R, such that Γ R n (U ) → Res(U ) as n → +∞ for any U ∈ Ω R , where the convergence is in the sense of the Attouch-Wets metric.That is, Henceforth, the radius R > 1 will remain fixed until stated otherwise.In order to prove Theorem 2.1, we start by studying the DtN maps.The goal is to recast the formula (1.2) into a form that can be implemented numerically.We introduce the orthonormal basis on ∂B R e n (θ which will be used frequently throughout the paper. 2.1.The outer map.The outer DtN map acts on a function φ ∈ H 1 (∂B R ) as In the orthonormal basis {e n } n∈Z the map M outer (k) has the explicit representation where H (1) ν denote the Hankel functions of the first kind.Using well-known identities for Bessel functions (cf.[21]), this can be rewritten as Moreover, it can be seen that ∼ z 2|n| as |n| → +∞, so that the matrix elements of M outer grow linearly in |n|.

2.2.
The inner map.We first show that the inner DtN map can be decomposed into a part that does not depend on U and a bounded part: Lemma 2.2.Let M inner,0 denote the free inner DtN map, i.e. the one with U = ∅.Then there exists a meromorphic operator-valued function Proof.For the inner DtN map, we first isolate the contribution of the resonator using the following calculation.
That is, S(k)φ is the harmonic extension of φ into B R .The regularity properties of S(k) follow from [24,Th. 4.21].Next choose any smooth radial function ρ satisfying (e.g.ρ(r) can be chosen as a piecewise polynomial in r; then the values ρ(r) can be computed from r in finitely many algebraic steps) and introduce the function In operator terms, this means that where H D denotes the Laplacian on L 2 (B R \ U ) with homogeneous Dirichlet boundary condition on ∂(B R \ U ) and where we have denoted Substituting the new representation for v into u = ρS(k)φ + v, we immediately obtain appearing in the last term in (2.6).The harmonic extension operator S(k) extends to a bounded operator [24,Th. 6.12]).Moreover, we have that T ρ : as a bounded operator.Finally, by the trace theorem, the normal derivative operator ∂ ν is continuous from is bounded.This completes the proof.
The decomposition in Lemma 2.2 allows us to reduce the study of M inner to that of M inner,0 .Next, we note that M inner,0 can in fact be represented as a diagonal operator in the basis {e n } n∈Z in a very similar manner to M outer .Indeed, one has where J ν denote the Bessel functions of the first kind.It can be seen that ∼ z 2|n| as |n| → +∞, so that the matrix elements of M inner,0 grow linearly in |n|.
2.3.Approximation procedure.From equations (2.1), (2.2) and (2.7) we have where , n ∈ Z Note that H, J , K are bounded and analytic in k.In order to determine whether ker(M inner + M outer ) = {0} we want to transform M inner + M outer into an operator of the form I + (compact).To this end, we introduce the bijective operator N := N + P 0 = diag(..., 3, 2, 1, 1, 1, 2, 3, ...), where P 0 denotes the projection onto the zeroth component.This new notation brings M inner + M outer into the form Because N . Therefore, its perturbation determinant det p is defined, where • denotes the ceiling function (cf.[11, Sec.XI.9]).Our task is reduced to finding the zeros k ∈ C − of the analytic function In order to approximate this determinant by something computable, we first prove that a square truncation of the operator matrix will converge in Schatten norm.
. .e n } be the orthogonal projection.Then there exists a constant C > 0 depending only on the set U such that Proof.The proof is given by a direct calculation: Our problem is therefore reduced to computing the matrix elements of the truncated operator This is done by performing a finite element procedure on B R \ U .This is the objective of the next section.
Remark 2.4.We remind the reader that the algorithm is assumed to have access to the values of Bessel and Hankel functions (recall (1.3)) so that the values of P n N − 1 2 (H + J )N − 1 2 P n − P 0 can be computed in finitely arithmetic operations.Hence we only need to compute P n N − 1 2 KN − 1 2 P n .

Matrix elements of the inner DtN map
The goal of this section is to compute the matrix elements where {e α } α∈Z is the orthonormal basis of L 2 (∂B R ) defined in the previous section and where we recall that . This is done by converting the integral on the boundary ∂B R to an integral on the interior, using Green's theorem.To this end, we extend the basis functions to the interior, by defining (in polar coordinates) where ρ was defined in (2.4), and denote Next, note that the integrand in the third term on the right hand side is explicitly given.Indeed, one can easily see that where J α are the Bessel functions of the first kind.Hence, the last integral in (3.2) can be explicitly approximated by quadrature, with a known bound for the error.In order to compute the first and second terms in the right hand side of (3.2), we need to approximate the function v α , which is the solution of the boundary value problem Note again that the right hand side of this equation is explicitly given in terms of (3.3).
3.1.Finite element approximation of v α .For notational convenience we write In the current subsection we prove: Proposition 3.1.Let v α be the solution of (3.4).For small discretization parameter h > 0 there exists a piecewise linear function v h α which is computable from a finite subset of (1.3) in finitely many algebraic steps, which satisfies the error estimate where C is independent of h and α, and where f α is defined in (3.1).
The reader should bear in mind that all constants obtained in our estimates will actually depend on k.However, these constants are uniformly bounded, as long as k 2 varies in a compact set separated from the spectrum of H D (recall that H D is the Laplacian on L 2 (O) with homogeneous Dirichlet boundary conditions on ∂O).Note that, since σ(H D ) is strictly positive and discrete, any compact subset of C − has this property.
Let h > 0 be our discretization parameter.We start by defining what we mean by the neighbors of a grid point j ∈ hZ 2 .Definition 3.2.Let j = (j 1 , j 2 ) ∈ hZ 2 .The set of neighbors of j is the set of all grid points contained in the cell [j 1 − h, j 1 + 2h] × [j 2 − h, j 2 + 2h] (cf. Figure 2).
Remark 3.3.The apparent asymmetry in the definition above is merely due to the fact that a point j defines a cell that is "northeast" of it, i.e. [j 1 , j

Define a uniform rectangular grid by
and note that L h can be completely constructed from the knowledge of χ U (j), j ∈ hZ 2 ∩ B R in finitely many steps.Next, define the open polygonal domain O h generated by L h as follows.
Note that as soon as h < 1 max(curv(∂O)) = 1 max(curv(∂U )) one will have O h ⊂ O and for all x ∈ ∂O h (cf. Figure 2) Figure 2. Illustration of eq.(3.5).The point j is included in L h , if and only if all the solid black points are in O.These solid black points are the neighbors of j.
Finally, we choose the obvious triangulation T h for O h obtained by splitting each square [0, h] 2 + j along its diagonal and denote by W h the corresponding P1 finite element space with zero boundary conditions, i.e.
Next, we will construct a finite element approximation for the solution of problem (3.4) starting from Céa's classical lemma: Then there exists C > 0 s.t.
for all h > 0.
We apply this lemma to our present situation by choosing Remark 3.5.Note that the finite element approximation v h α is not computable from the information in (1.3) in finitely many algebraic operations.An additional step of numerical integration is necessary to achieve this.We will revisit this issue at the end of the section.

Lemma 3.4 tells us that in order to estimate the approximation error
. For the latter, we have for any w h ∈ W h .To estimate the right hand side, we further decompose O h into a boundary layer and an interior part.More precisely, we write where O bdry h is the union of all cells having nonempty intersection with ∂O h , and Hence, for any w h ∈ W h we have Next, estimate the all three terms on the r.h.s. of (3.8) individually.To estimate the two last terms we will choose a particular "test function" u h ∈ W h .We first recall a classical interpolation estimate for linear finite elements [9, Th. 3.1.6]: Theorem 3.6 (Classical interpolation estimate).Let K be an element of T h and let Π h K denote the linear interpolation operator on H 2 (K), i.e.Π h K w is the affine function with (Π h K w)(j) = w(j) for all corners j of K. Then there exists C > 0 (independent of h, w, K) such that Now, we choose an explicit function u h ∈ W h as follows.
• For each finite element for boundary nodes i, u h (j) = v α (j), for interior nodes j.
Here a boundary node is any i ∈ L h ∩ ∂O h .Note that testing whether a node is a boundary node can be done in finitely many steps, because boundary nodes are precisely those nodes whose neighbors are not all contained in O.This defines a continuous, piecewise linear function in W h .We have the following error estimates.Lemma 3.7 (third term in (3.8)).One has the interior error estimate Proof.This follows immediately from Theorem 3.6.Lemma 3.8 (first term in (3.8)).One has the boundary layer error estimate Proof.We begin with v α − u h L 2 (O\O int h ) .Let x ∈ O \ O int h and note that similarly to eq. (3.5) there exists y ∈ ∂O with |x − y| ≤ 4 √ 2h.By [8, Ch. 9, Remark 7] this implies where the last line follows from the Gagliardo-Sobolev-Nirenberg inequality.Moreover, we note that by the definition of u h , one has . These facts allow us to estimate where we have used (3.10) in the last line.
Turning to the gradient, a similar estimate is obtained as follows.By definition of u h , one has In particular, for every K ⊂ O bdry h one has (where (3.12) was used in the fourth line and (3.10) in the last line).Hence the gradient estimate becomes (for h small enough).Finally, by elliptic regularity, we obtain the estimate where we recall that f α was defined in (3.1).Remark 3.9 (Numerical integration).As we noted in Remark 3.5, the approximate solution v h α is not computable from a finite subset of (1.3) in finitely many steps.However, the classical theory of finite elements shows that numerical integration does not change the error estimate (3.14) (see e.g.[9,Sec. 4.1]).This concludes the proof of Proposition 3.1.

3.2.
Matrix elements, revisited.Going back to (3.2) (and taking into account (3.3)), in order to compute the matrix elements K αβ , we have to approximate the integrals This will be a simple consequence of Prop.3.1 and Th.3.6.Let Π h denote the linear interpolation operator w.r.t. the triangulation T h , i.e.Π h u(j) = u(j) for all j ∈ L h and Π h u| K is an affine function for all elements K. Proposition 3.10.We have Proof.We begin with (3.16) and assume w.l.o.g.h < 1. Constants independent of h, α, β will be denoted C and their value may change from line to line.By the triangle and Hölder inequalities, we have where we have used Prop.3.1 & Th.3.6 in the second line, and Prop.3.1 & elliptic regularity in the third line.The estimate E β H 2 ≤ Cβ 2 used in the last line is evident from the definition of E β .We omit the calculations leading to (3.17) and (3.18), since they are entirely analogous.Remark 3.11.We note that (by a straightforward calculation in polar coordinates) the integral (3.15) vanishes unless α = β.In fact, one has where δ αβ denotes the Kronecker symbol.This one-dimensional integral can easily be approximated by a Riemann sum.While this fact is irrelevant to our theoretical treatment, it may improve performance and precision of numerical implementations.
Let us now go back to (3.2).Prop 3.10 allows us to define an approximate n × n matrix K h by whose matrix elements satisfy the error estimate The norms on the right hand side can be estimated by explicit calculation.One has Plugging these bounds into (3.20),we finally arrive at We summarize these results in the following Theorem 3.12.For any n ∈ N one has the operator norm error estimate * where C(k) is independent of n, h and bounded for k in a compact set disjoint from the real line.
Proof.This follows from the generalized Young inequality: (3.21), this immediately implies the assertion.Remark 3.13.We emphasize that all terms on the right hand side of (3.19) can be computed from the data provided in (1.3) in finitely many steps.Indeed, the mass and stiffness matrices m, s of the P1 finite element method can be computed exactly and satisfy for all φ, ψ ∈ W h .4. The algorithm 4.1.Proof of Theorem 2.1.From Theorem 3.12 (by choosing h ∼ n −9 ) and Lemma 2.3, we obtain a computable approximation K h(n) such that Using the Lipschitz continuity of perturbation determinants (cf.[27, Th. 6.5]) we conclude Then there exists C > 0, which is independent of k for k in a compact subset of C − , such that p .* We are abusing notation here, by identifying can only have finitely many zeros in Q and all are of finite order.That is, near a zero k 0 one has for some C, ν > 0 and for all k in a sufficiently small neighborhood of k 0 .Next, we define the R and Q dependent algorithm by The proof consists of two steps.First we prove that any convergent sequence k n ∈ Γ Q,R n (U ) necessarily converges to a zero of det p (I + A R ), and second we prove that for every zero k 0 of det p (I + A R ) there exists a sequence k n ∈ Γ Q,R n (U ) converging to k 0 .Together, these two facts imply Hausdorff convergence.
Let us begin with step 1. Assume that Hence, by Lemma 4.1 we have that as n → +∞.Hence we have Conversely, suppose that det(I + A R (k 0 )) = 0 for some k 0 ∈ Q.Then there exists a sequence Hence, by (4.2) we have n ν and employing Lemma 4.1 again, we obtain Clearly, for large enough n, the right hand side of (4.3) will be less than 1 log(n) , hence k n ∈ Γ Q,R n (U ) for n large enough.This completes the proof.
It remains to extend our argument from a single compact set Q to the entire complex plane.This is done via a diagonal-type argument.We choose a tiling of C − , where we start with a rectangle and then add rectangles in a counterclockwise spiral manner as shown in Figure 4.Note again that each individual rectangle Q i is well separated from {k | k 2 ∈ σ(H D )} ⊂ R. Next, we define our algorithm as follows.We let

4.2.
Proof of Theorem 1.9.It remains to extend the above algorithm from Ω R to Ω.In this section we will define a new algorithm, Γ n , based on the algorithm Γ R n constructed above.The main difficulty consists in "locating U " by testing only finitely many points.
The algorithm below uses two radii, R and r.R is used for running Γ R n , while r is successively increased to search for hidden components of U .The following pseudocode gives the definition of Γ n .n = 0: Let U ∈ Ω and initialise R = r = 1.n > 0: Consider the lattice L n = n −9 Z 2 ∩ B r (recall the choice h ∼ n −9 before eq.(4.1)).For every j ∈ L n , test whether j ∈ U .For all j ∈ L n ∩ U , compute |j|.
• Else, increment r by 1, set R := r and repeat the current step.This process generates increasing sequences, that we denote by R n and r n , and defines an algorithm Γ n : Ω → cl(C).Note that r n diverges to +∞, because it gets incremented by at least 1 in every step.
Lemma 4.3.The sequence {R n } n∈N is eventually constant and if Proof.The fact that {R n } n∈N is eventually constant follows immediately from the boundedness of U and the above pseudocode.Now let N ∈ N be as in the assertion.Assume for contradiction that U B R N −1 .Then, there exists x ∈ U with |x| > R N − 1, and, since But then, as soon as n > N is large enough such that n −9 < ε and r n > |x|, there would exist j ∈ L n with j ∈ B ε (x).According to the pseudocode, then, R n would be incremented by 1, contradicting the fact that R n = R N for all n > N .

By definition, the output of the algorithm Γ
where the second and third lines follow from Lemma 4.3 and the last line follows from the proof of Theorem 2.1.This completes the proof of Theorem 1.9.

Numerical Results
Although the algorithm described in Sections 3 and 4 was never designed for computational efficiency, we show in this section that it is possible to obtain surprisingly good numerical results with a MATLAB implementation which differs only in two minor respects: (a) The decomposition M inner = M inner,0 +K obtained in Lemma 2.2 introduces artificial singularities at the spectrum of H D (the Dirichlet Laplacian introduced in eq.(2.5)), which lead to numerical instabilities.For this reason, rather than approximating the solution v α of (3.4), we directly implemented a finite element approximation of the solution u of (2.3).The matrix elements of M inner can then be calculated via Green's formula, similarly to (3.2).(b) For reasons of performance, instead of using the triangulation procedure outlined in Section 3.1, we use the meshing tool Distmesh, cf.[26], to triangulate our domain.

5.1.
A circular resonator chamber.The first geometry we consider is that of a circular resonator chamber, connected to open space by a narrow channel.Figure 6 shows a logarithmic contour plot of the computed determinant | det(I + A n (k))| in the complex plane for n = 20.For the sake of comparison we included the Dirichlet eigenvalues of a closed resonator chamber in the image (red dots).One can see that next to each of the Dirichlet eigenvalues there is zero of the determinant, as we would expect.The resonance effects can in fact be seen on the level of the associated PDE. Figure 7 shows the finite element approximation of the solution u of (2.3) with right hand side φ = e α , once for k far away from a resonance, and once for k near a resonance.As can be seen, when k is near the resonance the wave penetrates into the chamber, whereas when k is far from a resonance the solution appears to be nearly 0 inside.[29,22,13].)This section contains a short study of the latter question, together with some numerical experiments.The approach outlined in Sections 2, 3 and 4 can be adapted to Neumann boundary conditions on the boundary of the obstacle U , given a set of stronger assumptions.In this section, we will briefly outline the ideas and present some numerical results.
Let now Ω N denote the class of bounded domains U B R−1 with C 2 boundary, on which we consider the Neumann Laplacians In addition, let us assume that each obstacle U comes with an associated regular parametrization We implemented the algorithm thus obtained (with the caveats (a), (b) mentioned at the beginning of the section) for an obstacle consisting of four circular holes, as shown in Figure 9.
This geometry had previously been studied in the context of water waves interacting with a circular array of cylinders [13].Figure 10 shows the output of our algorithm for four circular holes of radius 0.6, situated at the corners of a square of edge length 2.

Definition 1 . 1 .
Let us use the same symbols M inner (•), M outer (•) to denote the meromorphic continuations of the DtN maps to all of C. A number k ∈ C− := {z ∈ C | Im(z) < 0} in the lower half plane is called a resonance of H, if ker(M inner (k) + M outer (k)) = {0}.(1.2)Remark 1.2.Resonances in the sense of Definition 1.1 are well-defined.Indeed, M outer (k) is welldefined for all k ∈ C and the only poles of M inner (k) are equal to the Dirichlet eigenvalues of B R \ U .But these all lie on the real axis, which is excluded in Definition 1.1.

Lemma 3 . 4 (
Céa's Lemma, [9, Th. 2.4.1]).Let W h ⊂ H 1 0 (O) be a family of subspaces and assume the bilinear form a : H 1 0 (O) × H 1 0 (O) → C and linear functional f ∈ H −1 (O) satisfy the hypotheses of the Lax-Milgram theorem.Let v be the solution of the problem a(v, w) = f, w for all w ∈ H 1 0 (O) and v h be the solution of a(v h , w h ) = f, w h for all w h ∈ W h .(3.6)

Figure 3 Figure 3 .
Figure 3. Left: Sketch of O, O h and O bdry h .Right: Sketch of the chosen triangulation.

Figure 4 .
Figure 4. Tiling of the lower half plane

Figure 5 Figure 5 .
Figure 5. Example triangulation of a circular Helmholtz resonator used in our implementation.

Figure 6 .
Figure 6.Logarithmic contour plot of | det(I + An(k))| for n = 20 and opening width d = 1.3 in a strip below the real axis.Red dots: Dirichlet eigenvalues of the closed resonator chamber.

Figure 8 .
Figure 8. Top: Logarithmic contour plot of | det(I + An(k))| for n = 20 and opening width d = 1.0 in a strip below the real axis; Bottom left: zoomed in around first resonance; Bottom right: zoomed in around second resonance.Red dots: Dirichlet eigenvalues of the closed resonator chamber.

Figure 9 .
Figure 9. Triangulation of BR \ U for U consisting of four circular Neumann holes.

R 2 (
parametrized by its arc length) of its boundary.DefineΛ N := Λ ∪ U → γ U (x) x ∈ n U k=1 [a U k , b U k ] ,where Λ was defined in(1.3).Under this additional hypotheses, Sections 2, 3 and 4 (with modifications) show thatSCI(Ω N , Λ N , Res(•), (cl(C), d AW )) = 1.Indeed, the discussions in Sections 2 and 4 remain true with trivial changes for Neumann conditions on ∂U .The only major difference is in Section 3.1, where the finite element approximation of v α is constructed.A version of Prop.3.1 in the Neumann case can be obtained from our new assumptions as follows.(i) For each U , choose a uniform discretization of the intervals [a U k , b U k ].This induces an oriented discretization of ∂U via γ U .(ii) Using well known meshing algorithms (e.g.[7, Ch. 5]), construct a fitted mesh of B R \ U based on the boundary discretization from (i). (iii) Apply [2, Th. 4.1] to obtain the desired FEM error estimate.

Figure 10 .
Figure 10.Logarithmic contour plot of | det(I + An(k))| for n = 30 in a strip below the real axis

Table 1
summarizes the values of the first three eigenvalues and associated resonances in both cases (that is, with openings d = 1.0 and d = 1.3).One can observe that when the opening is narrower, the resonances are closer to their eigenvalue counterparts:

Table 1 .
Approximate numerical values of resonances and eigenvalues corresponding to the contour plots in Figures 6 and 8 respectively.5.2.Four circular Neumann holes.In many applications it is necessary to consider not onlyDirichlet, but also Neumann boundary conditions.An important real life situation in which this is the case is presented by resonance effects of water waves caused by submerged objects.A concrete question which has entailed a vast body of scientific literature is this: do the pillars of offshore structures, such as oil drilling platforms, introduce a positive interference of water waves that could damage the structure itself?At what frequencies are such positive interferences expected to occur?(cf.

Table 2 .
Approximate numerical values of resonances corresponding to the contour plot in Figure10.