Wavenumber-explicit analysis for the Helmholtz h-BEM: error estimates and iteration counts for the Dirichlet problem

Abstract

We consider solving the exterior Dirichlet problem for the Helmholtz equation with the h-version of the boundary element method (BEM) using the standard second-kind combined-field integral equations. We prove a new, sharp bound on how the number of GMRES iterations must grow with the wavenumber k to have the error in the iterative solution bounded independently of k as $k\rightarrow \infty$ when the boundary of the obstacle is analytic and has strictly positive curvature. To our knowledge, this result is the first-ever sharp bound on how the number of GMRES iterations depends on the wavenumber for an integral equation used to solve a scattering problem. We also prove new bounds on how h must decrease with k to maintain k-independent quasi-optimality of the Galerkin solutions as $k \rightarrow \infty$ when the obstacle is nontrapping.

1 Introduction

This paper is concerned with the wavenumber-explicit numerical analysis of boundary integral equations (BIEs) for the Helmholtz equation

$$\begin{aligned} \Delta u+k^2 u=0, \end{aligned}$$

(1.1)

where $k>0$ is the wavenumber, posed in the exterior of a 2- or 3-dimensional bounded obstacle ${\Omega }$ with Dirichlet boundary conditions on $\Gamma :=\partial {\Omega }$.

We consider the standard second-kind combined-field integral equation formulations of this problem: the so-called “direct” formulation (arising from Green’s integral representation)

$$\begin{aligned} A'_{k,\eta }v = f_{k,\eta } \end{aligned}$$

(1.2)

and the so-called “indirect” formulation (arising from an ansatz of layer potentials not related to Green’s integral representation)

$$\begin{aligned} A_{k,\eta }\phi = g_k, \end{aligned}$$

(1.3)

where

$$\begin{aligned} A'_{k,\eta }:= \frac{1}{2}I + D'_k -\mathrm{i}\eta S_k, \qquad A_{k,\eta }:= \frac{1}{2}I + D_k -\mathrm{i}\eta S_k, \end{aligned}$$

(1.4)

$\eta \in {\mathbb {R}}\setminus \{0\}$ is an arbitrary coupling parameter, $S_k$ is the single-layer operator, $D_k$ is the double-layer operator, and $D'_k$ is the adjoint double-layer operator (1.7), (1.8).

For simplicity of exposition, we focus on the direct Eq. (1.2), but the main results also hold for the indirect Eq. (1.3) (see Remark 1.20 below). The contribution to Eq. (1.2) from the Dirichlet boundary conditions is contained in the right-hand side $f_{k,\eta }$; our results are independent of the particular form of $f_{k,\eta }$, and so we can simplify the presentation by restricting attention to the particular exterior Dirichlet problem corresponding to scattering by a point source or plane wave, i.e. the sound-soft scattering problem (Definition 1.7 below).

We consider solving the Eq. (1.2) in ${L^2({\partial {\Omega }})}$ using the Galerkin method; this method seeks an approximation $v_N$ to the solution v from a finite-dimensional approximation space ${{{\mathcal {V}}}}_N$ (where N is the dimension, i.e. the total number of degrees of freedom). In the majority of the paper $\partial \Omega$ is $C^2$, in which case ${{{\mathcal {V}}}}_N$ will be the space of piecewise polynomials of degree p, for some fixed $p\ge 0$, on shape-regular meshes of diameter h, with h decreasing to zero; this is the so-called h–version of the Galerkin method, and we denote ${{{\mathcal {V}}}}_N$ and $v_N$ by ${{{\mathcal {V}}}}_h$ and $v_h$, respectively, and note that $N\sim h^{-(d-1)}$, where d is the dimension. To find the Galerkin solution $v_h$, one must solve a linear system of dimension N; in practice this is usually done using Krylov-subspace iterative methods such as the generalized minimal residual method (GMRES).

For the numerical analysis of this situation when k is large, there are now, roughly speaking, two main questions:

• Q1. How must h decrease with k in order to maintain accuracy of the Galerkin solution as $k\rightarrow \infty$?

• Q2. How does the number of GMRES iterations required to achieve a prescribed accuracy grow with k?

The goal of this paper is to prove rigorous results about these two questions, and then compare them with the results of numerical experiments.

We now give short summaries of the main results. These results depend on the choice of the coupling parameter $\eta$; for the results on Q1 we need $|\eta |\sim k$ and for the results on Q2 we need $\eta \sim k$, where we use the notation $a\sim b$ to mean that there exists $C_1, C_2>0$, independent of h and k, such that $C_1 b\le a \le C_2 b$. We also use the notation $a\lesssim b$ to mean that there exists $C>0$, independent of h and k, such that $a\le Cb$.

Summary of main results regarding Q1 and their context

Numerical experiments indicate that, in many cases, the condition $hk\lesssim 1$ is sufficient for the Galerkin method to be quasi-optimal (with the constant of quasi-optimality independent of k; i.e., (1.14) below holds); see [44, Section 5]. This feature can be described by saying that the h-BEM does not suffer from the pollution effect (in constrast to the h-FEM; see, e.g., [7, 52, Chapter 4]). The best existing result in the literature is that k-independent quasi-optimality of the Galerkin method applied to the integral Eq. (1.2) holds when $hk^{(d+1)/2}\lesssim 1$ for 2- and 3-d $C^{2,\alpha }$ obstacles that are star-shaped with respect to a ball [44, Theorem 1.4]. In this paper we improve this result by showing that the k-independent quasioptimality holds for 2-d nontrapping obstacles when $hk^{3/2}\lesssim 1$, for 3-d nontrapping obstacles when $hk^{3/2}\log k \lesssim 1$, and for 2- and 3-d smooth (i.e. $C^\infty$) convex obstacles with strictly positive curvature when $hk^{4/3}\lesssim 1$ (see Theorem 1.10 below).

The ideas behind the proofs of these results are summarised in Remark 1.13 below, but we highlight here that all the integral-operator bounds used in these arguments are sharp up to a factor of $\log k$. Therefore, to lower these thresholds on h for which quasi-optimality is proved, one would need to use different arguments than in the present paper. We also highlight that recent experiments by Baydoun and Marburg [10, 11, 60, 61] give examples of Helmholtz problems where the h-BEM suffers from a pollution effect, and therefore determining the sharp threshold on h for k-independent quasi-optimality to hold in general is an exciting open question.

Summary of main results regarding Q2 and their context There has been a large amount of research effort expended on understanding empirically how iteration counts for integral-equation formulations of scattering problems involving the Helmholtz or Maxwell equations depend on k; see, e.g, [1, 4, 15, 16, 81], and the references therein.

To our knowledge, however, there are no sharp k-explicit bounds in the literature, for any integral-equation formulation of a Helmholtz or Maxwell scattering problem, on the number of iterations GMRES requires to achieve a prescribed accuracy. The main reason, in this current setting of the Helmholtz exterior Dirichlet problem, is that the operator $A'_{k,\eta }$ is non-normal for all obstacles other than the circle and sphere [13, 14]. Therefore, for sufficiently-accurate discretisations, the Galerkin matrix of $A'_{k,\eta }$ is also non-normal, and one cannot use the well-known bounds on GMRES iterations in terms of the condition number (see, e.g., the review in [71, Section 6]).

In this paper, we prove that, for 2- and 3-d analytic obstacles with strictly positive curvature, the number of GMRES iterations growing like $k^{1/3}$ is sufficient to have the error in the iterative solution bounded independently of k (see Theorem 1.16 below). Numerical experiments in Sect. 5 show that the numbers of GMRES iterations for the sphere and an ellipsoid grow slightly less than $k^{1/3}$ ($k^{0.29}$ for the sphere and $k^{0.28}$ for an ellipsoid), and thus our bound is effectively sharp.

The ideas behind the proof are summarised in Remark 1.18 below. The focus of this paper is in proving results for the operator $A'_{k,\eta }$, i.e. the operator in the standard second-kind integral formulation, but we highlight in Remark 4.5 below how a bound on the number of GMRES iterations of $k^{1/2}$ when $d=2$ and $k^{1/2}\log k$ when $d=3$ can be obtained for a modification of $A'_{k,\eta }$, the so-called star-combined integral equation introduced in [74]. Moreover, whereas our bound on the number of iterations of $k^{1/3}$ for $A'_{k,\eta }$ holds for analytic obstacles with strictly positive curvature, the bounds for the star-combined operator hold for a much wider class of obstacles, namely piecewise-smooth Lipschitz obstacles that are star-shaped with respect to a ball.

Discussion of these results in the context of using semiclassical analysis in the numerical analysis of the Helmholtz equation

In the last 10 years, there has been growing interest in using results about the k-explicit analysis of the Helmholtz equation from semiclassical analysis (a branch of microlocal analysis) to design and analyse numerical methods for the Helmholtz equation.¹ The activity has occurred in, broadly speaking, four different directions:

• 1. The use of the results of Melrose and Taylor [64]–on the rigorous $k\rightarrow \infty$ asymptotics of the solution of the Helmholtz equation in the exterior of a smooth convex obstacle with strictly positive curvature—to design and analyse k-dependent approximation spaces for integral-equation formulations [5, 29,30,31,32, 41],

• 2. The use of the results of Melrose and Taylor [64], along with the work of Ikawa [53] on scattering from several convex obstacles, to analyse algorithms for multiple scattering problems [2, 33].

• 3. The use of bounds on the Helmholtz solution operator (also known as resolvent estimates) due to Vainberg [80] (using the propagation of singularities results of Melrose and Sjöstrand [63]) and Morawetz [66] to prove bounds on both $\Vert (A'_{k,\eta })^{-1}\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}$ and the inf-sup constant of the domain-based variational formulation [9, 22, 24, 73], and also to analyse preconditioning strategies [40].

• 4. The use of identities originally due to Morawetz [66] to prove coercivity of $A'_{k,\eta }$ [75] and to introduce new coercive formulations of Helmholtz problems [28, 42, 43, 65, 74].

The results of the present paper arise from a fifth direction, namely using estimates on the restriction of quasimodes of the Laplacian to hypersurfaces from [17, 26, 48, 77,78,79] to prove sharp k-explicit bounds on $S_k, D_k$ and $D'_k$ as operators from ${L^2({\partial {\Omega }})}$ to $H^1({\partial {\Omega }})$. We state these sharp k-explicit bounds in Sect. 2 below, and they are proved in the companion paper [39]. In the present paper, we use these new results, in conjunction with the results in Points 3 and 4 above, to obtain answers to Q1 and Q2.

1.1 Formulation of the problem

1.1.1 Geometric definitions

Let ${\Omega }\subset {\mathbb {R}}^d,$$d=2$ or 3,  be a bounded Lipschitz open set, such that the open complement ${\Omega _+}:= {\mathbb {R}}^d \setminus \overline{{\Omega }}$ is connected. Let $H^1_{\text {loc}}(\Omega _+)$ denote the set of functions v such that $\chi v \in H^1(\Omega _+)$ for every $\chi \in C^\infty _{\mathrm{comp}}(\overline{\Omega _+}):=\{ \chi |_{\Omega _+} : \,\chi \in C^\infty ({\mathbb {R}}^d) \text { is compactly supported}\}$. Let $\gamma ^{+}$ denote the trace operator from $\Omega _{+}$ to ${\partial {\Omega }}$. Let ${\mathbf {n}}$ be the outward-pointing unit normal vector to ${\Omega }$ (i.e. ${\mathbf {n}}$ points out of ${\Omega }$ and in to ${\Omega _+}$), and let $\partial _n^+$ denote the normal derivative trace operator from ${\Omega _+}$ to ${\partial {\Omega }}$ that satisfies $\partial _n^+ u= {\mathbf {n}}\cdot \gamma ^+(\nabla u)$ when $u\in H^2_{\text {loc}}({\Omega _+})$. (We also call $\gamma ^+ u$ the Dirichlet trace of u and $\partial _n^+ u$ the Neumann trace.)

Definition 1.1

(Star-shaped, and star-shaped with respect to a ball)

• (i) ${\Omega }$ is star-shaped with respect to the point${\mathbf {x}}_0\in {\Omega }$ if, whenever ${\mathbf {x}}\in {\Omega }$, the segment $[{\mathbf {x}}_0,{\mathbf {x}}]\subset {\Omega }$.

• (ii) ${\Omega }$ is star-shaped with respect to the ball$B_{a}({\mathbf {x}}_0)$ if it is star-shaped with respect to every point in $B_{a}({\mathbf {x}}_0)$.

• (iii) ${\Omega }$ is star-shaped with respect to a ball if there exists $a>0$ and ${\mathbf {x}}_0\in {\Omega }$ such that ${\Omega }$ is star-shaped with respect to the ball $B_{a}({\mathbf {x}}_0)$.

Definition 1.2

(Nontrapping) We say that ${\Omega }\subset {\mathbb {R}}^d, \,d=2, 3$ is nontrapping if ${\partial {\Omega }}$ is smooth ($C^\infty$) and, given R such that $\overline{{\Omega }}\subset B_R({\mathbf {0}})$, there exists a $T(R)<\infty$ such that all the billiard trajectories (in the sense of Melrose–Sjöstrand [63, Definition 7.20]) that start in ${\Omega _+}\cap B_R({\mathbf {0}})$ at time zero leave ${\Omega _+}\cap B_R({\mathbf {0}})$ by time T(R).

Definition 1.3

(Smooth hypersurface) We say that $\Gamma \subset {\mathbb {R}}^d$ is a smooth hypersurface if there exists $\widetilde{\Gamma }$ a compact embedded smooth $d-1$ dimensional submanifold of ${\mathbb {R}}^d$, possibly with boundary, such that $\Gamma$ is an open subset of $\widetilde{\Gamma }$, with $\Gamma$ strictly away from $\partial \widetilde{\Gamma }$, and the boundary of $\Gamma$ can be written as a disjoint union

$$\begin{aligned} \partial \Gamma =\left( \bigcup _{\ell =1}^n Y_\ell \right) \cup \Sigma , \end{aligned}$$

where each $Y_\ell$ is an open, relatively compact, smooth embedded manifold of dimension $d-2$ in $\widetilde{\Gamma }$, $\Gamma$ lies locally on one side of $Y_\ell$, and $\Sigma$ is closed set with $d-2$ measure 0 and $\Sigma \subset \overline{\bigcup _{l=1}^nY_l}$. We then refer to the manifold $\widetilde{\Gamma }$ as an extension of $\Gamma$.

For example, when $d=3$, the interior of a 2-d polygon is a smooth hypersurface, with $Y_i$ the edges and $\Sigma$ the set of corner points.

Definition 1.4

(Curved) We say a smooth hypersurface is curved if there is a choice of normal so that the second fundamental form of the hypersurface is everywhere positive definite.

Recall that the principal curvatures are the eigenvalues of the matrix of the second fundamental form in an orthonormal basis of the tangent space, and thus “curved” is equivalent to the principal curvatures being everywhere strictly positive (or everywhere strictly negative, depending on the choice of the normal).

Definition 1.5

(Piecewise smooth) We say that a hypersurface $\Gamma$ is piecewise smooth if $\Gamma =\cup _{i=1}^N \overline{\Gamma }_i$ where $\Gamma _i$ are smooth hypersurfaces and $\Gamma _i\cap \Gamma _j=\emptyset .$

Definition 1.6

(Piecewise curved) We say that a piecewise-smooth hypersurface $\Gamma$ is piecewise curved if $\Gamma$ is as in Definition 1.5 and each $\Gamma _j$ is curved.

1.1.2 The boundary value problem and integral equation formulation

Definition 1.7

(Sound-soft scattering problem) Given $k>0$ and an incident plane wave $u^I({\mathbf {x}}) = \exp (\mathrm{i}k{\mathbf {x}}\cdot {\mathbf {a}})$ for some ${\mathbf {a}}\in {\mathbb {R}}^d$ with $|{\mathbf {a}}|= 1$, find $u^S \in C^2(\Omega _+) \cap H^1_{\text {{loc}}}(\Omega _+)$ such that the total field $u:= u^I + u^S$ satisfies the Helmholtz equation (1.1) in ${\Omega _+}$, $\gamma ^+ u= 0$ on ${\partial {\Omega }}$, and $u^S$ satisfies the Sommerfeld radiation condition

$$\begin{aligned} \frac{\partial u^S}{\partial r}({\mathbf {x}}) - \mathrm{i}k \,u^S({\mathbf {x}}) = o\bigg (\frac{1}{r^{(d-1)/2}}\bigg ) \end{aligned}$$

as $r := |{\mathbf {x}}| \rightarrow \infty$, uniformly in ${\mathbf {x}}/r$.

The incident field in the sound-soft scattering problem of Definition 1.7 is a plane wave, but this could be replaced by a point source or, more generally, a solution of the Helmholtz equation in a neighbourhood of $\Omega$; see [19, Definition 2.11].

Obtaining the direct integral equation (1.2). If u satisfies the sound-soft scattering problem of Definition 1.7 then Green’s integral representation implies that

$$\begin{aligned} u({\mathbf {x}})= u^I({\mathbf {x}}) - \int _{\partial {\Omega }}\Phi _k({\mathbf {x}},{\mathbf {y}})\partial _n^+ u({\mathbf {y}})\,\mathrm{d}s({\mathbf {y}}), \quad {\mathbf {x}}\in {\Omega _+}, \end{aligned}$$

(1.5)

(see, e.g., [19, Theorem 2.43]), where $\Phi _k({\mathbf {x}},{\mathbf {y}})$ is the fundamental solution of the Helmholtz equation given by

$$\begin{aligned} \Phi _k({\mathbf {x}},{\mathbf {y}})=\displaystyle \frac{\mathrm{i}}{4}H_0^{(1)}\big (k|{\mathbf {x}}-{\mathbf {y}}|\big ), \,\,d=2,\quad \quad \Phi _k({\mathbf {x}},{\mathbf {y}}) = \frac{\mathrm{e}^{\mathrm{i}k |{\mathbf {x}}-{\mathbf {y}}|}}{4\pi |{\mathbf {x}}-{\mathbf {y}}|}, \,\,d=3 \end{aligned}$$

(note that we have chosen the sign of $\Phi _k({\mathbf {x}},{\mathbf {y}})$ so that $-(\Delta +k^2)\Phi _k({\mathbf {x}},{\mathbf {y}})= \delta ({\mathbf {x}}-{\mathbf {y}})$). Taking the exterior Dirichlet and Neumann traces of (1.5) on ${\partial {\Omega }}$ and using the jump relations for the single- and double-layer potentials (see, e.g., [19, Equation 2.41]) we obtain the integral equations

$$\begin{aligned} S_k \partial _n^+ u= \gamma ^+ u^I \quad \text { and }\quad \left( \frac{1}{2}I + D_k'\right) \partial _n^+ u= \partial _n^+ u^I, \end{aligned}$$

(1.6)

where $S_k$ and $D'_k$ are the single- and adjoint-double-layer operators defined by

$$\begin{aligned} S_k \phi ({\mathbf {x}}) :=&\int _{\partial {\Omega }}\Phi _k({\mathbf {x}},{\mathbf {y}}) \phi ({\mathbf {y}})\,\mathrm{d}s({\mathbf {y}}), \quad D'_k \phi ({\mathbf {x}}) := \int _{\partial {\Omega }}\frac{\partial \Phi _k({\mathbf {x}},{\mathbf {y}})}{\partial n({\mathbf {x}})} \phi ({\mathbf {y}})\,\mathrm{d}s({\mathbf {y}}), \end{aligned}$$

(1.7)

for $\phi \in {L^2({\partial {\Omega }})}$ and $x\in {\partial {\Omega }}$. Later we will also need the definition of the double-layer potential,

$$\begin{aligned} D_k \phi ({\mathbf {x}}) := \int _{\partial {\Omega }}\frac{\partial \Phi _k({\mathbf {x}},{\mathbf {y}})}{\partial n({\mathbf {y}})} \phi ({\mathbf {y}})\,\mathrm{d}s({\mathbf {y}}) \quad \text { for }\phi \in {L^2({\partial {\Omega }})}\,\text { and }\, {\mathbf {x}}\in {\partial {\Omega }}. \end{aligned}$$

(1.8)

The first equation in (1.6) is not uniquely solvable when $-k^2$ is a Dirichlet eigenvalue of the Laplacian in ${\Omega }$, and the second equation in (1.6) is not uniquely solvable when $-k^2$ is a Neumann eigenvalue of the Laplacian in ${\Omega }$ (see, e.g., [19, Theorem 2.25]). The standard way to resolve this difficulty is to take a linear combination of the two equations, which yields the integral Eq. (1.2) where $A'_{k,\eta }$ is defined by (1.4),

$$\begin{aligned} f_{k,\eta }:= \partial _n^+ u^I- \mathrm{i}\eta \, \gamma ^+ u^I, \end{aligned}$$

(1.9)

and we use the notation that $v:=\partial _n^+ u$ (this makes denoting the Galerkin solution below easier, since we then have $v_h$ instead of $(\partial _n^+ u)_h$).

The space ${L^2({\partial {\Omega }})}$ is a natural space for the practical solution of second-kind integral equations since it is self-dual, and, for $\eta \in {\mathbb {R}}\setminus \{0\}$, $A'_{k,\eta }$ is a bounded invertible operator from ${L^2({\partial {\Omega }})}$ to itself [19, Theorem 2.27]. Furthermore the right-hand side $f_{k,\eta }$ is in ${L^2({\partial {\Omega }})}$ (since $u^I\in C^\infty (\overline{{\Omega _+}})$) and thus we consider the Eq. (1.2) as an equation in ${L^2({\partial {\Omega }})}$.

The Galerkin method Given a finite-dimensional approximation space ${{{\mathcal {V}}}}_N \subset {L^2({\partial {\Omega }})}$, the Galerkin method for the integral Eq. (1.2) is

$$\begin{aligned} \text {find } v_N \in {{{\mathcal {V}}}}_N \text { such that }\, \big (A'_{k,\eta }v_N, w_N\big )_{{L^2({\partial {\Omega }})}}=\big (f_{k,\eta },w_N\big )_{{L^2({\partial {\Omega }})}} \, \text { for all } w_N\in {{{\mathcal {V}}}}_N.\nonumber \\ \end{aligned}$$

(1.10)

Let ${{{\mathcal {V}}}}_N= \mathop {\mathrm{span}}\{\phi _i: i=1,\ldots ,N\}$, let $v_N \in {{{\mathcal {V}}}}_N$ be equal to $\sum _{j=1}^N V_j \phi _j$, and define ${\mathbf {v}}\in {\mathbb {C}}^N$ by ${\mathbf {v}}:= (V_j)_{j=1}^N$. Then, with ${\mathbf {A}}_{ij}:= (A'_{k,\eta }\phi _j,\phi _i)_{{L^2({\partial {\Omega }})}}$ and ${\mathbf {f}}_i:= (f_{k,\eta },\phi _i)_{{L^2({\partial {\Omega }})}}$, the Galerkin method (1.10) is equivalent to solving the linear system ${\mathbf {A}}{\mathbf {v}}= {\mathbf {f}}$.

We consider the h–version of the Galerkin method, and we then denote ${{{\mathcal {V}}}}_N$ and $v_N$ by ${{{\mathcal {V}}}}_h$ and $v_h$ respectively. The main results for Q1 and Q2 will be stated under the following assumption on ${{{\mathcal {V}}}}_h$.

Assumption 1.8

(Assumptions on${{{\mathcal {V}}}}_h$) ${{{\mathcal {V}}}}_h$ is a space of piecewise polynomials of degree p for some fixed $p\ge 0$ on shape-regular meshes of diameter h, with h decreasing to zero (see, e.g., [70, Chapter 4] for specific realisations). Furthermore

• (a) if $w\in H^1({\partial {\Omega }})$ then

  $$\begin{aligned} \min _{w_h \in {{{\mathcal {V}}}}_h}\left\| w - w_h\right\| _{{L^2({\partial {\Omega }})}} \lesssim h\left\| w\right\| _{H^1({\partial {\Omega }})}, \end{aligned}$$

  (1.11)
• (b)
  $$\begin{aligned} \left\| w_h\right\| ^2_{{L^2({\partial {\Omega }})}} \sim h^{d-1} \left\| {\mathbf {w}}\right\| _2^2, \end{aligned}$$

  (1.12)

  where $\Vert \cdot \Vert _2$ denotes the $l_2$ (i.e. euclidean) vector norm.

Remark 1.9

(For what situations is Assumption 1.8 proved?) Part (a) is proved for subspaces consisting of piecewise-constant basis functions in [70, Theorem 4.3.19] when ${\partial {\Omega }}$ is a polyhedron or curved (in the sense of Assumptions 4.3.17 and 4.3.18, respectively, in [70]) and in [76, Theorem 10.4] when ${\Omega }$ is a piecewise-smooth Lipschitz domain. Part (a) is proved for subspaces consisting of continuous piecewise-polynomials of degree $p\ge 1$ (in the sense of [70, Definition 4.1.36]) in [70, Theorem 4.3.28].

Part (b) is proved for subspaces consisting of piecewise-linear basis function in [76, Lemma 10.5] when ${\partial {\Omega }}$ is piecewise-smooth and Lipschitz, and for more general subspaces in [70, Theorem 4.4.7].

1.2 Statement of the main results and discussion

The results concerning Q1 are stated in Sect. 1.2.1, and the results concerning Q2 are stated in Sect. 1.2.2.

1.2.1 Results concerning Q1

Theorem 1.10

(Sufficient conditions for the Galerkin method to be quasi-optimal) Let u be the solution of the sound-soft scattering problem of Definition 1.7 and let $v:=\partial _n^+ u$. Let $|\eta |\sim k$, and let ${{{\mathcal {V}}}}_h$ satisfy Part (a) of Assumption 1.8.

• (a) If either (i) ${\Omega }$ is nontrapping, or (ii) ${\Omega }$ is star-shaped with respect to a ball and ${\partial {\Omega }}$ is $C^{2,\alpha }$ and piecewise smooth, then given $k_0>0$, there exists a $C>0$ (independent of k and h) such that if

  $$\begin{aligned} h k^{3/2} \le C, \quad d=2,\quad h k^{3/2}\log k \le C, \quad d=3, \end{aligned}$$

  (1.13)

  then the Galerkin equations (1.10) have a unique solution which satisfies

  $$\begin{aligned} \left\| v-v_h\right\| _{{L^2({\partial {\Omega }})}} \lesssim \min _{w_h\in {{{\mathcal {V}}}}_h} \left\| v-w_h\right\| _{{L^2({\partial {\Omega }})}} \end{aligned}$$

  (1.14)

  for all $k\ge k_0$.

• (b) In case (ii) above, if additionally ${\partial {\Omega }}$ is piecewise curved, then given $k_0>0$, there exists a $C>0$ (independent of k and h) such that if

  $$\begin{aligned} h k^{4/3}\log k \le C, \quad d=2,3 \end{aligned}$$

  (1.15)

  then (1.14) holds.

• (c) If ${\Omega }$ is convex and ${\partial {\Omega }}$ is $C^\infty$ and curved then given $k_0>0$, there exists a $C>0$ (independent of k and h) such that if

  $$\begin{aligned} h k^{4/3} \le C, \quad d=2,3 \end{aligned}$$

  (1.16)

  then (1.14) holds.

Having established quasi-optimality, it is then natural to ask how the best approximation error $\min _{w_h\in {{{\mathcal {V}}}}_h} \left\| v-w_h\right\| _{{L^2({\partial {\Omega }})}}$ depends on k, h, and $\Vert v\Vert _{L^2({\partial {\Omega }})}$.

Theorem 1.11

(Bounds on the best approximation error) Let u be the solution of the sound-soft scattering problem of Definition 1.7 and let $v:=\partial _n^+ u$. Let ${{{\mathcal {V}}}}_h$ satisfy Assumption 1.8.

• (a) If ${\partial {\Omega }}$ is $C^{2,\alpha }$ and piecewise smooth, then, given $k_0>0$,

  $$\begin{aligned} \min _{w_h \in {{{\mathcal {V}}}}_h}\left\| v-w_h\right\| _{{L^2({\partial {\Omega }})}} \lesssim h A(k) \left\| v\right\| _{{L^2({\partial {\Omega }})}} \end{aligned}$$

  (1.17)

  with $A(k)= k^{5/4}\log k$, for all $k\ge k_0$.

• (b) If ${\partial {\Omega }}$ is piecewise curved, then, given $k_0>0$, (1.17) holds with $A(k)= k^{7/6}\log k$, for all $k\ge k_0$.

• (c) If ${\Omega }$ is convex and ${\partial {\Omega }}$ is $C^\infty$ and curved, then, given $k_0>0$, (1.17) holds with $A(k)= k$, for all $k\ge k_0$.

Combining Theorems 1.10 and 1.11 we can obtain bounds on the relative error of the Galerkin method. For brevity, we only state the ones corresponding to cases (a) and (c) in Theorems 1.10 and 1.11.

Corollary 1.12

(Bound on the relative errors in the Galerkin method) Let u be the solution to the sound-soft scattering problem, let $|\eta |\sim k$, and let ${{{\mathcal {V}}}}_h$ satisfy Part (a) of Assumption 1.8.

• (a) If either (i) ${\Omega }$ is nontrapping, or (ii) ${\Omega }$ is star-shaped with respect to a ball and ${\partial {\Omega }}$ is $C^{2,\alpha }$ and piecewise smooth, then given $k_0>0$, there exists a $C>0$ (independent of k and h) such that if h and k satisfy (1.13) then the Galerkin equations (1.10) have a unique solution which satisfies

  $$\begin{aligned} \frac{\left\| v-v_h\right\| _{{L^2({\partial {\Omega }})}}}{\left\| v\right\| _{{L^2({\partial {\Omega }})}}} \lesssim {\left\{ \begin{array}{ll} k^{-1/4}\log k ,&{} d=2,\\ k^{-1/4}, &{} d=3, \end{array}\right. } \end{aligned}$$

  for all $k\ge k_0$.

• (b) If ${\Omega }$ is convex and ${\partial {\Omega }}$ is $C^\infty$ and curved, then given $k_0>0$, there exists a $C>0$ (independent of k and h) such that if $h k^{4/3} \le C$ the Galerkin equations (1.10) have a unique solution which satisfies

  $$\begin{aligned} \frac{\left\| v-v_h\right\| _{{L^2({\partial {\Omega }})}}}{\left\| v\right\| _{{L^2({\partial {\Omega }})}}} \lesssim \frac{1}{k^{1/3}} \end{aligned}$$

  for all $k\ge k_0$.

Remark 1.13

(The main ideas behind the proofs of Theorems 1.10 and 1.11 ) The proof of Theorem 1.10 uses the classic projection-method analysis of second-kind integral equations (see, e.g., [6, Chapter 3]), with $A'_{k,\eta }$ treated as a compact perturbation of $\frac{1}{2}I$. In [44], this argument was used to reduce the question of finding k-explicit bounds on the mesh threshold h for k-independent quasi-optimality to hold to finding k-explicit bounds on

$$\begin{aligned} \Vert S_k\Vert _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}, \quad \Vert D'_k\Vert _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}, \quad \text { and }\quad \Vert (A'_{k,\eta })^{-1}\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}. \end{aligned}$$

We use the new, sharp bounds on the first two of these norms from [39], quoted here as Theorem 2.1, and the sharp bounds on the third of these norms from [9, Theorem 1.13] (for nontrapping obstacles) and [22, Theorem 4.3] (for obstacles that are star-shaped with respect to a ball).

The bounds of Theorem 1.11 are proved by showing that

$$\begin{aligned} \left\| v\right\| _{H^1({\partial {\Omega }})}\lesssim A(k) \left\| v\right\| _{{L^2({\partial {\Omega }})}}, \end{aligned}$$

(1.18)

and then using the approximation theory result (1.11). The bound (1.18) is obtained from the integral equation (1.2) using the second-kind-structure of the equation and the ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$ bounds on $S_k$ and $D'_k$ from Theorem 2.1.

Remark 1.14

(Comparison to previous results) Theorems 1.10 and 1.11 and Corollary 1.12 sharpen previous results in [44]: the mesh thresholds for quasi-optimality in Theorem 1.10 are sharper than the corresponding ones in [44], and the results are valid for a wider class of obstacles.

This sharpening is due to the new, sharp bounds on ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$ norms of $S_k$, $D_k$, and $D'_k$ from [39], and the widening of the class of obstacles is due to the bound on $\Vert (A'_{k,\eta })^{-1}\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}$ for nontrapping obstacles from [9, Theorem 1.13]. In more detail: Theorem 1.4 of [44] is the analogue of our Theorem 1.10 except that the former is only valid when ${\Omega }$ is star-shaped with respect to a ball and $C^{2,\alpha }$ and the mesh threshold is $hk^{(d+1)/2}\le C$. Comparing this result to Theorem 1.10 we see that we’ve sharpened the threshold in the $d=3$ case, expanded the class of obstacles to nontrapping ones, and added the additional results (b) and (c). Theorem 1.11 on the best approximation error is again proved using the ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$-bounds from [39] and thus we see similar improvements over the corresponding theorem in [44] ( [44, Theorem 1.3]).

As discussed in Remark 1.13, both the present paper and [44] use the classic projection-method argument to obtain k-explicit results about quasi-optimality of the h-BEM. There are two other sets of results about quasi-optimality of the h-BEM in the literature:

• (a) results that use coercivity [29, 74, 75], and

• (b) results that give sufficient conditions for quasi-optimality to hold in terms of how well the spaces ${{{\mathcal {V}}}}_h$ approximate the solution of certain adjoint problems [8, 57, 62].

These two sets of results are discussed in detail in [44, pp. 181–182] and [44, Section 4.2] respectively, and neither give results as strong as those in Theorem 1.10.

Finally, in this paper we have only considered the h-BEM; a thorough k-explicit analysis of the hp-BEM for the exterior Dirichlet problem was conducted in [57] and [62]. In particular, this analysis, combined with the bound on $\Vert (A'_{k,\eta })^{-1}\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}$ for nontrapping obstacles from [9, Theorem 1.13], proves that k-independent quasi-optimality can be obtained for nontrapping obstacles through a choice of h and p that keeps the total number of degrees of freedom proportional to $k^{d-1}$ [57, Corollaries 3.18 and 3.19].

Remark 1.15

(How sharp are the quasioptimality results?) Numerical experiments in [44, Section 5] show that for a wide variety of obstacles (including certain mildly-trapping obstacles) the h-BEM is quasi-optimal with constant independent of k (i.e. (1.14) holds), when $hk\sim 1$. The closest we can get to proving this is the result for strictly convex obstacles in Theorem 1.10 part (c), with the threshold being $hk^{4/3}\le C$. The recent results of Baydoun and Marburg [10, 11, 60, 61], however, give examples of cases where $hk\sim 1$ is not sufficient to keep the error bounded as $k\rightarrow \infty$.

1.2.2 Result concerning Q2

We now consider solving the linear system ${\mathbf {A}}{\mathbf {v}}= {\mathbf {f}}$ with the generalised minimum residual method (GMRES) introduced by Saad and Schultz in [69]; for details of the implementation of this algorithm, see, e.g., [45, 68].

Theorem 1.16

(A bound on the number of GMRES iterations) Let ${\Omega }$ be a 2- or 3-d convex obstacle whose boundary ${\partial {\Omega }}$ is analytic and curved. Let ${{{\mathcal {V}}}}_h$ satisfy Part (b) of Assumption 1.8, let the Galerkin matrix corresponding to (1.10) be denoted by ${\mathbf {A}}$, and consider GMRES applied to ${\mathbf {A}}{\mathbf {v}}= {\mathbf {f}}$

There exist constants $\eta _0>0$ and $k_0>0$ (with $\eta _0=1$ if $\Omega$ is a ball) such that if $k\ge k_0$ and $\eta _0 k \le \eta \lesssim k$, then, given $0<\varepsilon <1$, there exists a C (independent of k, $\eta$, and $\varepsilon$) such that if

$$\begin{aligned} m\ge C k^{1/3}\log \left( \frac{12}{\varepsilon }\right) , \end{aligned}$$

(1.19)

then the mth GMRES residual ${\mathbf {r}}_m:= {\mathbf {A}}{\mathbf {v}}_m -{\mathbf {f}}$ satisfies

$$\begin{aligned} \frac{\left\| {\mathbf {r}}_m\right\| _2}{\left\| {\mathbf {r}}_0\right\| _2}\le \varepsilon , \end{aligned}$$

where $\Vert \cdot \Vert _2$ denotes the $l_2$ (i.e. euclidean) vector norm.

In other words, Theorem 1.16 states that, for convex, analytic, curved $\Omega$, the number of iterations growing like $k^{1/3}$ is a sufficient condition for GMRES to maintain accuracy as $k\rightarrow \infty$.

Remark 1.17

(How sharp is the result of Theorem 1.16?) Numerical experiments in Sect. 5 show that for the sphere the number of GMRES iterations grows like $k^{0.29}$, and for an ellipsoid they grow like $k^{0.28}$. The bound in Theorem 1.16 is therefore effectively sharp (at least for the range of k considered in the experiments).

Remark 1.18

(The main ideas behind the proof of Theorem 1.16) The two ideas behind Theorem 1.16 are that:

• (a) A sufficient (but not necessary) condition for iterative methods to be well behaved is that the numerical range (also known as the field of values) of the matrix is bounded away from zero, and in this case the Elman estimate [34, 35] and its refinement due to Beckermann, Goreinov, and Tyrtyshnikov [12] can be used to bound the number of GMRES iterations in terms of (i) the distance of the numerical range to the origin, and (ii) the norm of the matrix.

• (b) When ${\Omega }$ is convex, $C^3$, piecewise analytic, and ${\partial {\Omega }}$ is curved, [75] proved that $A'_{k,\eta }$ is coercive for sufficiently large k (with $\eta \sim k$). The k-dependence of the coercivity constant, along with the k-dependence of $\Vert A'_{k,\eta }\Vert _{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}$ then give the information needed about the numerical range of the Galerkin matrix ${\mathbf {A}}$ required in (a).

Remark 1.19

(Comparison to previous results) The bound $m > rsim k^{2/3}$ when ${\partial {\Omega }}$ is a sphere was stated in [75, Section 1.3]; this bound was obtained using the original Elman estimate (see Remark 4.4 below), and the fact that the sharp bound $\Vert A'_{k,\eta }\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}\lesssim k^{1/3}$ was known for the circle and sphere; see [19, Section 5.4]. To our knowledge, there are no other k-explicit bounds in the literature on the number of GMRES iterations required to achieve a prescribed accuracy for a Helmholtz BIE. The closest related work is [23], which uses a second-kind integral equation to solve the Helmholtz equation in a half-plane with an impedance boundary condition. The special structure of this integral equation allows a two-grid iterative method to be used, and [23] proves that there exists $C>0$ such that if $k h\le C$, then, after seven iterations, the difference between the solution and the Galerkin solution computed via the iterative method is bounded independently of k and h.

Remark 1.20

(Translating the results to the indirect equation (1.3)) Instead of using Green’s integral representation (1.5) to formulate the sound-soft scattering problem as the integral equation (1.2), one can pose the ansatz that the scattered field satisfies

$$\begin{aligned} u^S({\mathbf {x}})=\int _{\partial {\Omega }}\frac{\partial \Phi _k({\mathbf {x}},{\mathbf {y}})}{\partial n({\mathbf {y}})} \phi ({\mathbf {y}})\,\mathrm{d}s({\mathbf {y}}) - \mathrm{i}\eta \int _{\partial {\Omega }}\Phi _k({\mathbf {x}},{\mathbf {y}})\phi ({\mathbf {y}})\,\mathrm{d}s({\mathbf {y}}) \end{aligned}$$

for ${\mathbf {x}}\in {\Omega _+}$, $\phi \in {L^2({\partial {\Omega }})}$, and $\eta \in {\mathbb {R}}\setminus \{0\}$. Imposing the boundary condition $\gamma ^+u^S=-\gamma ^+ u^I$ on ${\partial {\Omega }}$ and using the jump relations for the single- and double-layer potentials leads to the integral equation (1.3) where $A_{k,\eta }$ is defined by (1.4) and $g= - \gamma ^+ u^I$. One can show that $A_{k,\eta }$ and $A'_{k,\eta }$ are adjoint with respect to the real-valued ${L^2({\partial {\Omega }})}$ inner product (see, e.g., [19, Equation 2.37, Remark 2.24, §2.6]), and so their norms are equal, the norms of their inverses are equal, and if one is coercive then so is the other (with the same coercivity constant). These facts imply that the results of Theorems 1.10 and Theorem 1.16 hold for the indirect equation (1.3).

The bounds on the best approximation error in Theorem 1.11 hold for the indirect equation (1.3) with (a) $A(k) = k^{3/2}$ for $d=2$, $A(k)=k^{3/2}\log k$ for $d=3$, (b) $A(k)=k^{4/3}\log k$, and (c) $A(k)=k^{4/3}$. These powers of k are all slightly higher than those for the direct equation; the reason for this is that we have more information about the unknown in the direct equation (since it is $\partial _n^+ u$) than about the unknown $\phi$ in the indirect equation. Indeed, one can express $\phi$ in terms of the difference of solutions to interior and exterior boundary value problems–see [19, p. 132]–but it is harder to make use of this fact than for the direct equation.

Remark 1.21

(Translating the results to the general exterior Dirichlet problem) The results of Theorems 1.10 and 1.16 are independent of the right-hand side of the integral equation (1.2), and therefore hold for the general Dirichlet problem with Dirichlet data in $H^1({\partial {\Omega }})$ (this assumption is needed so that $A'_{k,\eta }$ can still be considered as an operator on ${L^2({\partial {\Omega }})}$; see, e.g., [19, Section 2.6]). The results of Theorem 1.11 and Corollary 1.12, however, do not immediately hold for the general Dirichlet problem, since they use the particular form of the right-hand side in (1.9).

Outline of the paper In Sect. 2 we recap the sharp ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$ bounds from the companion paper [39]. In Sect. 3 we prove Theorems 1.10 and 1.11 (the results concerning Q1). In Sect. 4 we prove Theorem 1.16 (the result concerning Q2), and then in Sect. 5 we give numerical experiments showing that Theorem 1.16 is effectively sharp in its k-dependence.

2 Recap of the ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$ bounds from [39]

The following result was proved in [39, Theorem 2.10]. In stating this result, we use the weighted $H^1({\partial {\Omega }})$ norm

$$\begin{aligned} \left\| w\right\| _{H^1_k({\partial {\Omega }})}^2:= k^{-2}\left\| \nabla _{{\partial {\Omega }}}w\right\| _{{L^2({\partial {\Omega }})}}^2+ \left\| w\right\| _{{L^2({\partial {\Omega }})}}^2, \end{aligned}$$

(2.1)

in contrast to the usual $H^1({\partial {\Omega }})$ norm

$$\begin{aligned} \left\| w\right\| _{H^1({\partial {\Omega }})}^2:=\left\| \nabla _{{\partial {\Omega }}}w\right\| _{{L^2({\partial {\Omega }})}}^2+ \left\| w\right\| _{{L^2({\partial {\Omega }})}}^2, \end{aligned}$$

where $\nabla _{{\partial {\Omega }}}$ is the surface gradient operator on ${\partial {\Omega }}$ (see, e.g., [19, p. 276]); we note that the use of such weighted norms is standard in both the semiclassical and numerical analysis of the Helmholtz equation, and reflects the fact that we expect to incur a power of k every time we take a derivative of a Helmholtz solution; see, e.g., [65, Remark 3.8].

Theorem 2.1

(Bounds on $\Vert S_k\Vert _{{L^2({\partial {\Omega }})}\rightarrow H^1_k({\partial {\Omega }})}$, $\Vert D_k\Vert _{{L^2({\partial {\Omega }})}\rightarrow H^1_k({\partial {\Omega }})}$, $\Vert D'_k\Vert _{{L^2({\partial {\Omega }})}\rightarrow H^1_k({\partial {\Omega }})}$)

Let $\Omega$ be a bounded Lipschitz open set such that the open complement ${\Omega _+}:= {\mathbb {R}}^d\setminus \overline{\Omega }$ is connected.

• (a) If ${\partial {\Omega }}$ is a piecewise-smooth hypersurface (in the sense of Definition 1.5), then, given $k_0>1$, there exists $C>0$ (independent of k) such that

  $$\begin{aligned} \Vert S_k\Vert _{L^2({\partial {\Omega }})\rightarrow H_{k}^1({\partial {\Omega }})} \le C\,k^{-1/2}\,\log k. \end{aligned}$$

  (2.2)

  for all $k\ge k_0$. Moreover, if ${\partial {\Omega }}$ is piecewise curved (in the sense of Definition 1.6), then, given $k_0>1$, the following stronger estimate holds for all $k\ge k_0$

  $$\begin{aligned} \left\| S_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow {H}_k^1({\partial {\Omega }})} \le C k^{-2/3}\log k. \end{aligned}$$

  (2.3)

• (b) If ${\partial {\Omega }}$ is a piecewise smooth, $C^{2,\alpha }$ hypersurface, for some $\alpha >0$, then, given $k_0>1$, there exists $C>0$ (independent of k) such that

  $$\begin{aligned} \left\| D_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow {H}_k^1({\partial {\Omega }})} + \left\| D'_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow {H}_k^1({\partial {\Omega }})} \le C k^{1/4}\log k. \end{aligned}$$

  Moreover, if ${\partial {\Omega }}$ is piecewise curved, then, given $k_0>1$, there exists $C>0$ (independent of k) such that the following stronger estimates hold for all $k\ge k_0$

  $$\begin{aligned} \left\| D_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1_k({\partial {\Omega }})} + \left\| D'_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1_k({\partial {\Omega }})} \lesssim k^{1/6}\log k. \end{aligned}$$

• (c) If ${\Omega }$ is convex and ${\partial {\Omega }}$ is $C^\infty$ and curved (in the sense of Definition 1.4) then, given $k_0>1$, there exists C such that, for $k\ge k_0$,

  $$\begin{aligned}&\displaystyle \left\| S_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H_k^1({\partial {\Omega }})} \le C k^{-2/3}, \\&\left\| D_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H_k^1({\partial {\Omega }})}+\left\| D'_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H_k^1({\partial {\Omega }})} \le C. \end{aligned}$$

The requirement in Part (b) of Theorem 2.1 that ${\partial {\Omega }}$ is $C^{2,\alpha }$ arises since this is the regularity required of ${\partial {\Omega }}$ for $D_k$ and $D'_k$ to map ${L^2({\partial {\Omega }})}$ to $H^1({\partial {\Omega }})$; see [54, Theorem 4.2], [27, Theorem 3.6].

The bounds in Theorem 2.1 contain k-explicit ${L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}$ bounds on $S_k, D_k$ and $D'_k$. These ${L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}$ bounds were originally proved in [47, Appendix A] and [36] (and the realisation that these ${L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}$ bounds could be extended to ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$ bounds was the motivation for [39]).

Remark 2.2

(Sharpness of the bounds in Theorem 2.1) In [39, Section 3] it is shown that, modulo the factor $\log k$, all of the bounds in Theorem 2.1 are sharp (i.e. the powers of k in the bounds are optimal). The sharpness (modulo the factor $\log k$) of the ${L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}$ bounds contained in Theorem 2.1 was proved in [47, Section  A.2–A.3]. Earlier work in [18, Section 4] proved the sharpness of some of the ${L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}$ bounds in 2-d; we highlight that [39, Section 3] and [47, Section  A.2–A.3] contain the appropriate generalisations to multidimensions of some of the arguments of [18, Section 4] (in particular [18, Theorems 4.2 and 4.4]).

Remark 2.3

(Sharp bounds on $S_k$ when $d=2$) When $d=2$ and ${\partial {\Omega }}$ is Lipschitz, the sharp bound

$$\begin{aligned} \Vert S_k\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}\lesssim k^{-1/2} \end{aligned}$$

(2.4)

was proved using the Riesz–Thorin interpolation theorem in [18, Theorem 3.3] and by the Schur test in [38, Theorem 6]. Similarly, the sharp bound

$$\begin{aligned} \left\| S_k\right\| _{{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}} \lesssim k^{1/2} \end{aligned}$$

(2.5)

was proved using the Riesz–Thorin interpolation theorem in [44, Theorem 1.6].

3 Proofs of Theorems 1.10, 1.11 (the results concerning Q1)

3.1 Proof of Theorem 1.10

The heart of the proof of Theorem 1.10 is the following lemma.

Lemma 3.1

There exists a $\tilde{C}>0$ such that under the condition

$$\begin{aligned} h \left\| D_k'-\mathrm{i}\eta S_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})} \Vert (A'_{k,\eta })^{-1}\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}} \le \widetilde{C} \end{aligned}$$

(3.1)

the Galerkin equations (1.10) have a unique solution satisfying (1.14).

The presence of $\Vert (A'_{k,\eta })^{-1}\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}$ in (3.1) means that before proving Theorem 1.10 using Lemma 3.1 we need to recall the following bounds on $\Vert (A'_{k,\eta })^{-1}\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}$.

Theorem 3.2

([22, Theorem 4.3], [9, Theorem 1.13]) If $|\eta |\sim k$ and either${\Omega }$ is star-shaped with respect to a ball and $C^2$ in a neighbourhood of almost every point on $\Gamma$or${\Omega }$ is nontrapping, then, given $k_0>0$, $\Vert (A'_{k,\eta })^{-1}\Vert _{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}\lesssim 1$ for all $k\ge k_0$.

Proof of Theorem 1.10 using Lemma 3.1

Using the triangle inequality, a sufficient condition for (3.1) to hold is

$$\begin{aligned}&h\left( \left\| D_k'\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})} + |\eta |\left\| S_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}\right) \nonumber \\&\quad \Vert (A'_{k,\eta })^{-1}\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}} \le \widetilde{C}. \end{aligned}$$

(3.2)

In [39, Remark 2.22] it is shown that the ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$ norms of $D_k'$ and $S_k$ are maximised in different regions of phase space, and thus we do not lose anything by using the triangle inequality, i.e., (3.2) is no less sharp than (3.1) in terms of k-dependence.

The mesh thresholds (1.13), (1.15), (1.16) then follow from using the bound $\Vert (A'_{k,\eta })^{-1}\Vert _{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}\lesssim 1$ from Theorem 3.2 and the different bounds on $\Vert D_k'\Vert _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}$ and $\Vert S_k\Vert _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}$ in Theorem 2.1 (recalling the definition of the $H^1_k({\partial {\Omega }})$ norm in (2.1)), apart from when $d=2$ when we use the bound on $S_k$ (2.5) instead of (2.2). $\square$

To prove Theorem 1.10 we therefore only need to prove Lemma 3.1. This was proved in [44, Corollary 4.1], but since the proof is short we repeat it here for completeness.

We first introduce some notation: let $P_h$ denote the orthogonal projection from ${L^2({\partial {\Omega }})}$ onto ${{{\mathcal {V}}}}_h$ (see, e.g, [6, Section 3.1.2]); then the Galerkin equations (1.10) are equivalent to the operator equation

$$\begin{aligned} P_h A'_{k,\eta }v_h = P_h f_{k,\eta }. \end{aligned}$$

(3.3)

The proof requires us to treat $A'_{k,\eta }$ as a (compact) perturbation of the identity, and thus we let $L_{k,\eta }:= D_k' - \mathrm{i}\eta S_k$. Furthermore, to make the notation more concise, we let $\lambda =1/2$. Therefore, the left-hand side of (3.3) becomes $(\lambda I + P_h L_{k,\eta })v_h$, and the question of existence of a solution to (3.3) boils down to the invertibility of $(\lambda I + P_h L_{k,\eta })$. Note also that, using the $P_h$ notation, the best approximation error on the right-hand side of (1.14) is $\Vert (I-P_h)v\Vert _{{L^2({\partial {\Omega }})}}$.

The heart of the proof of Lemma 3.1 is the following lemma.

Lemma 3.3

If

$$\begin{aligned} \left\| (I - P_h) L_{k,\eta }\right\| _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}} \Vert (A'_{k,\eta })^{-1}\Vert _{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}\le \frac{\delta }{1+ \delta } \end{aligned}$$

(3.4)

for some $\delta >0$, then the Galerkin equations have a unique solution, $v_h$, which satisfies the quasi-optimal error estimate

$$\begin{aligned} \left\| v - v_h\right\| _{{L^2({\partial {\Omega }})}} \le \lambda (1+ \delta ) \Vert (A'_{k,\eta })^{-1}\Vert _{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}\Vert (I-P_h)v\Vert _{{L^2({\partial {\Omega }})}}. \end{aligned}$$

(3.5)

Proof of Lemma 3.1 using Lemma 3.3

By the polynomial-approximation result (1.11),

$$\begin{aligned} \left\| (I-P_h)L_{k,\eta }\right\| _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}} \lesssim h \left\| L_{k,\eta }\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}. \end{aligned}$$

Therefore, choosing, say, $\delta =1$, we find that there exists a $\tilde{C}>0$ such that (3.1) implies that (3.4) holds. $\square$

Thus, to prove Theorem 1.10, we only need to prove Lemma 3.3.

Proof of Lemma 3.3

Since

$$\begin{aligned} \lambda I + P_h L_{k,\eta }&= \lambda I + L_{k,\eta } - (I-P_h)L_{k,\eta }\\&= \left( \lambda I + L_{k,\eta }\right) \left( I - \left( \lambda I+ L_{k,\eta }\right) ^{-1} (I-P_h)L_{k,\eta }\right) , \end{aligned}$$

if

$$\begin{aligned} \big \Vert \left( \lambda I+L_{k,\eta }\right) ^{-1} (I-P_h)L_{k,\eta }\big \Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}<1, \end{aligned}$$

then $(\lambda I + P_h L_{k,\eta })$ is invertible using the classical result that $I-A$ is invertible if $\left\| A\right\| <1$. In this abstract setting $\Vert (I-A)^{-1}\Vert \le (1-\left\| A\right\| )^{-1}$, and thus if (3.4) holds we have

$$\begin{aligned} \big \Vert (\lambda I + P_h L_{k,\eta })^{-1}\big \Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}&\le \big \Vert \left( \lambda I +L_{k,\eta }\right) ^{-1}\big \Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}} \frac{1}{1 - \delta /(1+\delta ) },\nonumber \\&= (1 + \delta ) \, \big \Vert (\lambda I + L_{k,\eta })^{-1}\big \Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}. \end{aligned}$$

(3.6)

Writing the direct equation as $(\lambda I + L_{k,\eta })v=f$ and the Galerkin equation as $(\lambda I + P_h L_{k,\eta })v_h=P_h f$, we have

$$\begin{aligned} v - v_h = v - (\lambda I + P_h L_{k,\eta })^{-1} P_h f&= (\lambda I + P_h L_{k,\eta })^{-1}(\lambda v - P_h(f - L_kv)) \nonumber \\ {}&= \lambda \left( \lambda I+ P_h L_{k,\eta }\right) ^{-1} (I- P_h ) v, . \end{aligned}$$

(3.7)

and the result (3.5) follows from using the bound (3.6) in (3.7). $\square$

Remark 3.4

(Is there a better choice of $\eta$ than $|\eta |\sim k$?) Theorem 1.10 is proved under the assumption that $|\eta |\sim k$. This choice of $\eta$ is widely recommended from studies of the condition number of $A'_{k,\eta }$; see [19, Chapter 5] for an overview of these. From (3.2) we see that the best choice of $\eta$, from the point of view of obtaining the least-restrictive threshold for k-independent quasi-optimality, will minimise the k-dependence of

$$\begin{aligned} \left( \left\| D_k'\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})} + |\eta |\left\| S_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}\right) \Vert (A'_{k,\eta })^{-1}\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}. \end{aligned}$$

There does not yet exist a rigorous proof that $|\eta |\sim k$ minimises this quantity, but [9, Section 7.1] outlines exactly the necessary results still to prove.

3.2 Proof of Theorem 1.11

Proof of Theorem 1.11

By the polynomial-approximation result (1.11), we only need to prove that the bound (1.18) hold with the different functions A(k). The idea is to take the $H^1$ norm of the integral equation (1.2) and then use the ${{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}$ and ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$ bounds contained in Theorems 2.1.

Taking the $H^1$ norm of (1.2) and using the notation that $A'_{k,\eta }= \frac{1}{2}I + L_{k,\eta }$ and $v:=\partial _n^+ u$, as in the proof of Theorem 1.10 above, we have that

$$\begin{aligned} \left\| v\right\| _{H^1({\partial {\Omega }})} \lesssim \left\| L_{k,\eta }\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}\left\| v\right\| _{{L^2({\partial {\Omega }})}} + \left\| f_{k,\eta }\right\| _{H^1({\partial {\Omega }})}. \end{aligned}$$

In this inequality, $\eta$ is just a parameter that appears in $L_{k,\eta }$ and $f_{k,\eta }$, with the equation holding for all values of $\eta$; in other words, the unknown $v(=\partial _n^+ u)$ does not depend on the value of $\eta$. We now seek to minimise the k-dependence of $\Vert L_{k,\eta }\Vert _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}$. Looking at the k-dependence of the ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$-bounds on $S_k$ and $D'_k$ in Theorem 2.1, we see that, under each of the different geometric set-ups, the best choice is $\eta =0$, and thus

$$\begin{aligned} \left\| v\right\| _{H^1({\partial {\Omega }})} \lesssim \left\| D'_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}\left\| v\right\| _{{L^2({\partial {\Omega }})}} + k^2 \end{aligned}$$

(3.8)

where we have explicitly worked out the k-dependence of $\Vert f_{k,\eta }\Vert _{H^1({\partial {\Omega }})}$ using the definition (1.9).

Taking the $L^2$ norm of (1.2) (with $\eta =0$), and noting that $\Vert f_{k,\eta }\Vert _{{L^2({\partial {\Omega }})}}\sim k$, we have that

$$\begin{aligned} \big (1 + \left\| D'_k\right\| _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}\big ) \left\| v\right\| _{{L^2({\partial {\Omega }})}} > rsim k. \end{aligned}$$

(3.9)

Using (3.9) in (3.8), we have

$$\begin{aligned} \left\| v\right\| _{H^1({\partial {\Omega }})} \lesssim \left( \left\| D'_k\right\| _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}+ k\big (1 + \left\| D'_k\right\| _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}\big )\right) \left\| v\right\| _{{L^2({\partial {\Omega }})}}.\nonumber \\ \end{aligned}$$

(3.10)

Since the bounds on the ${L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})$-norm of $D'_k$ in Theorem 2.1 are one power of k higher that the ${L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}$-bounds in Theorem 2.1, using these norm bounds in (3.10) results in the bound $\Vert v\Vert _{H^1({\partial {\Omega }})}\lesssim A(k)\Vert v\Vert _{{L^2({\partial {\Omega }})}}$ with the functions of A(k) as in the statement of theorem (and equal to the right-hand sides of the bounds on $\Vert D_k'\Vert _{{L^2({\partial {\Omega }})}\rightarrow H^1({\partial {\Omega }})}$ in Theorem 2.1). $\square$

4 Proofs of Theorem 1.16 (the result concerning Q2)

To prove Theorem 1.16 we need to recall (i) the result about coercivity of $A'_{k,\eta }$ when ${\Omega }$ is convex, $C^3$, piecewise analytic, and curved from [75], and (ii) the refinement of the Elman estimate in [12].

Theorem 4.1

(Coercivity of $A'_{k,\eta }$ for ${\Omega }$ convex, $C^3$, piecewise analytic, and curved [75]) Let ${\Omega }$ be a convex domain in either 2- or 3-d whose boundary, ${\partial {\Omega }}$, is curved and is both $C^3$ and piecewise analytic. Then there exist constants $\eta _0>0$, $k_0>0$ (with $\eta _0=1$ when $\Omega$ is a ball) and a function of k, $\alpha _{k}>0$, such that for $k\ge k_0$ and $\eta \ge \eta _0 k$,

$$\begin{aligned} \big |\big (A'_{k,\eta }\phi , \phi \big )_{{L^2({\partial {\Omega }})}} \big |\ge \alpha _{k}\Vert \phi \Vert ^2_{{L^2({\partial {\Omega }})}}\quad \text { for all }\phi \in {L^2({\partial {\Omega }})}, \end{aligned}$$

(4.1)

where

$$\begin{aligned} \alpha _{k} = \frac{1}{2}- {{{\mathcal {O}}}}\big ( k^{-2/3}\log k\big ) \quad \text { as }k \rightarrow \infty . \end{aligned}$$

(4.2)

In stating this result we have used the bound on $S_k$ (2.3) and [75, Remark 3.3] to get the asymptotics (4.2). The fact that $\eta _0=1$ when ${\Omega }$ is a ball follows from [74, Corollary 4.8].

Theorem 4.2

(Refinement of the Elman estimate [12]) Let ${\mathbf {A}}$ be a matrix with $0\notin W({\mathbf {A}})$, where $W({\mathbf {A}}):= \big \{ ({\mathbf {A}}{\mathbf {v}}, {\mathbf {v}}) : {\mathbf {v}}\in {\mathbb {C}}^N, \Vert {\mathbf {v}}\Vert _2=1\big \}$ is the numerical range of ${\mathbf {A}}$. Let $\beta \in [0,\pi /2)$ be defined such that

$$\begin{aligned} \cos \beta = \frac{\mathrm {dist}\big (0, W({\mathbf {A}})\big )}{\Vert {\mathbf {A}}\Vert _{2}}, \end{aligned}$$

and let $\gamma _\beta$ be defined by

$$\begin{aligned} \gamma _\beta := 2 \sin \left( \frac{\beta }{4-2\beta /\pi }\right) . \end{aligned}$$

(4.3)

Suppose the matrix equation ${\mathbf {A}}{\mathbf {v}}= {\mathbf {f}}$ is solved using GMRES, and let ${\mathbf {r}}_m:= {\mathbf {A}}{\mathbf {v}}_m - {\mathbf {f}}$ be the m-th GMRES residual. Then

$$\begin{aligned} \frac{\Vert {\mathbf {r}}_m\Vert _{2}}{\Vert {\mathbf {r}}_0\Vert _{2}} \le \left( 2 + \frac{2}{\sqrt{3}}\right) \big (2+ \gamma _\beta \big ) \,\gamma _\beta ^m. \end{aligned}$$

(4.4)

When we apply the estimate (4.4) to ${\mathbf {A}}$, we find that $\beta = \pi /2-\delta$, where $\delta = \delta (k)$ is such that $\delta \rightarrow 0$ as $k\rightarrow \infty$. We therefore specialise the result (4.4) to this particular situation in the following corollary.

Corollary 4.3

If $\beta = \pi /2-\delta$ with $0<\delta <\delta _0$, then there exists $C_1>0$ and $\delta _1>0$ (both independent of $\delta$) such that, for $0<\varepsilon <1$,

$$\begin{aligned} \text {if} \quad m\ge \frac{C_1}{\delta }\log \left( \frac{12}{\varepsilon }\right) \quad \text { then } \quad \frac{\Vert {\mathbf {r}}_m\Vert _{D}}{\Vert {\mathbf {r}}_0\Vert _{D}}\le \varepsilon \end{aligned}$$

(4.5)

for all $0<\delta <\delta _1$.

That is, choosing $m > rsim \delta ^{-1}$ is sufficient for GMRES to converge in an $\delta$-independent way as $\delta \rightarrow 0$.

Proof of Corollary 4.3

If $\beta = \pi /2-\delta$, with $\delta \rightarrow 0$, then $\cos \beta =\sin \delta = \delta + {{{\mathcal {O}}}}(\delta ^3)$ as $\delta \rightarrow 0$. From the definition of the convergence factor $\gamma _\beta$, (4.3), we have

$$\begin{aligned} \gamma _\beta:= & {} 2 \sin \left( \frac{\beta }{4-2\beta /\pi }\right) \nonumber \\= & {} 2 \sin \left( \frac{\pi }{6} -\frac{4\delta }{9}+{{{\mathcal {O}}}}(\delta ^2)\right) = 1 -\frac{4\delta }{3\sqrt{3}} + {{{\mathcal {O}}}}(\delta ^2)\quad \text { as }\,\, \delta \rightarrow 0, \end{aligned}$$

(4.6)

and then

$$\begin{aligned} \log \gamma _\beta = -\frac{4\delta }{3\sqrt{3}} + {{{\mathcal {O}}}}(\delta ^2) \quad \text { as }\,\, \delta \rightarrow 0, \end{aligned}$$

and so there exist $C_2>0$ and $\delta _1>0$ such that

$$\begin{aligned} \gamma _\beta ^m = \mathrm{e}^{m\log \gamma _\beta }\le \mathrm{e}^{-m\delta /C_2} \quad \text { for all }0< \delta \le \delta _1. \end{aligned}$$

The bound (4.5) then follows from (4.4) since $(2 + 2/\sqrt{3})(2+ \gamma _\beta )<3(2 + 2/\sqrt{3})<12$.

Remark 4.4

(Comparison of (4.4) with the original Elman estimate) The estimate

$$\begin{aligned} \frac{\Vert {\mathbf {r}}_m\Vert _{2}}{\Vert {\mathbf {r}}_0\Vert _{2}} \le \sin ^m \beta \end{aligned}$$

(4.7)

was essentially proved in [34, 35] (see also the review [71, Section 6] and the references therein). When $\beta = \pi /2-\delta$, the convergence factor in (4.7) is

$$\begin{aligned} \sin \beta = \cos \delta = 1- \frac{\delta ^2}{2}+{{{\mathcal {O}}}}(\delta ^4); \end{aligned}$$

by comparing this to (4.6) we can see that (4.7) is indeed a weaker bound.

Proof of Theorem 1.16

The set up of the Galerkin method in §1.1 implies that, for any $v_N, w_N\in {{{\mathcal {V}}}}_N$, $( A'_{k,\eta }v_N, w_N)_{{L^2({\partial {\Omega }})}} = ({\mathbf {A}}{\mathbf {v}},{\mathbf {w}})_2$, where $(\cdot ,\cdot )_2$ denotes the euclidean inner product on $l^2$. Therefore, the continuity of $A'_{k,\eta }$ and the norm equivalence (1.12) implies that

$$\begin{aligned} |({\mathbf {A}}{\mathbf {v}}, {\mathbf {w}})_2| \lesssim \left\| A'_{k,\eta }\right\| _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}} h^{d-1} \left\| {\mathbf {v}}\right\| _2 \left\| {\mathbf {w}}\right\| _2 \quad \text { for all }{\mathbf {v}}, {\mathbf {w}}\in {\mathbb {C}}^N. \end{aligned}$$

(4.8)

Furthermore, if $A'_{k,\eta }$ is coercive with coercivity constant $\alpha _{k,\eta }$, i.e., (4.1) holds, then

$$\begin{aligned} |({\mathbf {A}}{\mathbf {v}}, {\mathbf {v}})_2| > rsim \alpha _{k,\eta }h^{d-1} \left\| {\mathbf {v}}\right\| ^2_2 \quad \text { for all }{\mathbf {v}}\in {\mathbb {C}}^N. \end{aligned}$$

(4.9)

The bounds (4.8) and (4.9) together imply that the ratio $\cos \beta$ in (4.7) satisfies

$$\begin{aligned} \cos \beta > rsim \frac{\alpha _{k,\eta }}{\Vert A'_{k,\eta }\Vert _{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}. \end{aligned}$$

Since ${\Omega }$ is $C^\infty$ and curved, the bound $\Vert A'_{k,\eta }\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}\lesssim k^{1/3}$ follows from the bounds in Theorem 2.1 (recalling that $\eta _0 k\le \eta \lesssim k$). Since ${\partial {\Omega }}$ is piecewise analytic, $C^3$, and curved, from Theorem 4.1 there exists a $k_0>0$ such that $\alpha _{k,\eta }\sim 1$ for all $k\ge k_0$. Combining these two bounds we have $\cos \beta > rsim k^{-1/3}$ for all $k\ge k_0$ and thus Corollary 4.3 holds with $\delta \sim k^{-1/3}$ for all $k\ge k_0$; the result (1.19) then follows from (4.5).

Note that the assumption in the theorem that ${\partial {\Omega }}$ is analytic comes from the fact that if ${\partial {\Omega }}$ is both piecewise analytic and $C^\infty$, then ${\partial {\Omega }}$ must be analytic, where the notion of piecewise analyticity in Theorem 4.1 is inherited from [25, Definition 4.1]. $\square$

Remark 4.5

(The star-combined operator) The bound on the number of iterations in Theorem 1.16 crucially depends on the coercivity result of Theorem 4.1. Although numerical experiments in [14] indicate that $A'_{k,\eta }$ is coercive, uniformly in k, for a wider class of obstacles that those in Theorem 4.1, this has yet to be proved.²

Nevertheless, there does exist an integral operator that (i) can be used to solve the sound-soft scattering problem, and (ii) is provable coercive for a wide class of obstacles. Indeed, the star-combined operator${\mathscr {A}}_k$, introduced in [74] and defined by

$$\begin{aligned} {\mathscr {A}}_k:= \big ({\mathbf {x}}\cdot {\mathbf {n}}({\mathbf {x}})\big )\left( \frac{1}{2}I + D'_k\right) + {\mathbf {x}}\cdot \nabla _{\partial {\Omega }}S - \mathrm{i}\eta S_k \end{aligned}$$

(where $\nabla _{{\partial {\Omega }}}$ is the surface gradient operator on ${\partial {\Omega }}$; see, e.g., [19, p. 276]), has the following two properties:

• (i) if u solves the sound-soft scattering problem, then

  $$\begin{aligned} {\mathscr {A}}_k\partial _n^+ u= {\mathbf {x}}\cdot \gamma ^+(\nabla u^I)- \mathrm{i}\eta \gamma ^+ u^I \end{aligned}$$

  (4.10)

  [74, Lemma 4.1] (see also [19, Theorem 2.36]), and

• (ii) if ${\Omega }$ is a 2- or 3-d Lipschitz obstacle that is star-shaped with respect to a ball and $\eta := k|{\mathbf {x}}| + \mathrm{i}(d-1)/2$, then

  $$\begin{aligned} \mathfrak {R}\big ({\mathscr {A}}_k\phi ,\phi \big )_{{L^2({\partial {\Omega }})}} \ge \frac{1}{2}\mathop {\mathrm{ess} \inf }_{{\mathbf {x}}\in {\partial {\Omega }}}\big ({\mathbf {x}}\cdot {\mathbf {n}}({\mathbf {x}})\big )>0 \end{aligned}$$

  for all $k>0$ [74, Theorem 1.1]. The refinement of the Elman estimate in Theorem 4.2 can therefore be used to prove results about the number of iterations required when GMRES is applied to the Galerkin discretisation of (4.10). Since the coercivity constant of the star-combined operator is independent of k, the k-dependence of the analogue of the bound (1.19) for ${\mathscr {A}}_k$ rests on the bounds on $\Vert {\mathscr {A}}_k\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}$. For convex ${\Omega }$ with smooth and curved ${\partial {\Omega }}$, Theorem 2.1 implies that $\Vert {\mathscr {A}}_k\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}\lesssim k^{1/3}$, and we therefore obtain the same bound on m as for $A'_{k,\eta }$ (i.e. (1.19)). For general piecewise-smooth Lipschitz obstacles that are star-shaped with respect to a ball, Theorem 2.1 combined with the bounds (2.4) and (2.5) shows that $\Vert {\mathscr {A}}_k\Vert _{{{L^2({\partial {\Omega }})}\rightarrow {L^2({\partial {\Omega }})}}}\lesssim k^{1/2}$ when $d=2$ and $\lesssim k^{1/2}\log k$ when $d=3$. Corollary 4.3 then implies that $m > rsim k^{1/2}$ for $d=2$ and $m > rsim k^{1/2}\log k$ for $d=3$. Recall that GMRES always converges in at most N steps (in exact arithmetic), and when $h\sim 1/k$ we have that $N\sim k^{d-1}$; these bounds on m are therefore nontrivial.

5 Numerical experiments concerning Q2

The main purpose of this section is to show that the $k^{1/3}$ growth in the number of iterations given by Theorem 1.16 is effectively sharp.

Details of the scattering problems considered

We solve the sound-soft scattering problem of Definition 1.7 with ${\mathbf {a}}= (1,0,0)$ (i.e the incident plane wave propagates in the $x_1$-direction), using the direct integral equation (1.2) and the Galerkin method (1.10). The subspace ${{{\mathcal {V}}}}_h$ is taken to be piecewise constants on a shape regular mesh, and the meshwidth h is taken to be $2\pi /(10 k)$, i.e. we are choosing ten points per wavelength. We solve the resulting linear system with GMRES, with tolerance $1\times 10^{-5}$. We consider two obstacles:

• 1. ${\Omega }$ the unit sphere, and

• 2. ${\Omega }$ the ellipsoid with semi-principal axes of lengths 3, 1, and 1 (in the $x_1$-, $x_2$-, and $x_3$-directions respectively.

The computations were carried out using version 3.0.3 of the BEM++ library [72] on one node of the “Balena” cluster at the University of Bath. The cluster consists of Intel Xeon E5-2650 v2 (Ivybridge, 2.60 GHz) CPUs and the used node had 512GB of main memory. BEM++ was compiled with version 5.2 of the GNU C compiler and the Python code was run under Anaconda 2.3.0.

Numerical results

Tables 1 and 2 displays the number of degrees of freedom, number of iterations required for GMRES to converge, and time taken to converge, with $\eta =k$, and with ${\Omega }$ the sphere or ellipsoid. The difference between Tables 1 and 2 is that, in the first, k starts as 2 and then doubles until it equals 128, and in the second, k starts as 3 and then doubles until it equals 96; we performed the second set of experiments when the $k=128$ run for the ellipsoid failed to complete. Figure 1 plots the iteration counts from both tables and compares them to the $k^{1/3}$ rate from Theorem 1.16 (the graph is plotted on a $\log$-$\log$ scale so that a dependence $\#_{\text {iterations}}\sim k^\alpha$ appears as a straight line with gradient $\alpha$).

Table 1 With ${\Omega }$ the sphere or ellipsoid and $\eta =k$, the number of degrees of freedom, number of iterations required for GMRES to converge (with tolerance $1\times 10^{-5}$), and time taken to converge, when GMRES is applied to the Galerkin matrix corresponding to the direct integral equation (1.2), starting with $k=4$ and then doubling until $k=128$

Table 2 Same as Table 1 but for a different range of k

Fig. 1

The number of iterations required for GMRES to converge (with tolerance $1\times 10^{-5}$) when GMRES is applied to the Galerkin matrix corresponding to the direct integral equation (1.2) with $\eta =k$, and with ${\Omega }$ the sphere or ellipsoid, and the values of k from Tables 1 and 2 . The $k^{1/3}$ rate is the upper bound on the rate guaranteed by Theorem 1.16

We see from Fig. 1 that the $k^{1/3}$ growth predicted by Theorem 1.16 appears to be effectively sharp. Indeed, the plot of the iterations for the ellipsoid becomes roughly linear from $k=12$ onwards, and estimating the slope of this line using the numbers of iterations at $k=12$ and $k=96$ we have that the $\#_{\text {iterations}}\sim k^{0.28}$. Using the numbers of iterations at $k=12$ and $k=96$ to estimate the rate of growth for the sphere we have that $\#_{\text {iterations}}\sim k^{0.29}$.

Finally, Table 3 compares the iteration counts and times for the sphere when $\eta =k$ and when $\eta =-k$. We see that, for every value of k considered, the number of iterations when $\eta =-k$ is much greater than when $\eta =k$. Table 3 only goes up to $k=32$, since the $k=64$ run for the sphere with $\eta =-k$ did not complete.

Table 3 With ${\Omega }$ the sphere and $\eta =k$ or $\eta =-k$, the number of iterations required for GMRES to converge (with tolerance $1\times 10^{-5}$) and time taken to converge, when GMRES is applied to the Galerin matrix corresponding to the direct integral equation (1.2)

Remark 5.1

(The link between Table 3 and the recent work of Marburg [58, 59])

We performed the experiment in Table 3 because, in the engineering-acoustics literature, Marburg recently considered collocation discretisations of the direct integral equation for the Neumann problem (i.e. the Neumann-analogue of equation (1.2)) and showed that the analogue of the choice $\eta =k$ leads to much slower growth than the analogue of the choice $\eta =-k$ [58, 59].

A heuristic explanation for this dependence of the number of iterations on the sign of $\eta$ is essentially contained in the work of Levadoux and Michielsen [55, 56], and Antoine and Darbas [3]. In our setting of using the operator $A'_{k,\eta }$ to solve the exterior Dirichlet problem, the key points are that

• 1. the ideal $\mathrm{i}\eta$ should approximate the Dirichlet-to-Neumann (DtN) map in ${\Omega _+}$, and

• 2. $\mathrm{i}k$ is a better approximation to the DtN map than $-\mathrm{i}k$ (at least for smooth convex obstacles).

Regarding 1: taking the Dirichlet trace of Green’s integral representation (written with general Dirichlet data, not just data coming from a plane wave as in (1.5)), and using the jump relations for the single- and double-layer potentials (see, e.g., [19, Equation 2.41]) we find that

$$\begin{aligned} \gamma ^+ u= -S_k(\partial _n^+ u)+ \left( \frac{1}{2}I+ D_k\right) \gamma ^+ u. \end{aligned}$$

Rearranging this equation, and introducing the notation $P_{\mathrm{{DtN}}}^+$ for the exterior Dirichlet-to-Neumann map for solutions of the Helmholtz equation in ${\Omega _+}$ satisfying the Sommerfeld radiation condition, we find that

$$\begin{aligned} I = \frac{1}{2}I + D_k- S_k P_{\mathrm{{DtN}}}^+. \end{aligned}$$

(5.1)

Green’s second identity implies that, for $\phi ,\psi \in H^{1/2}({\partial {\Omega }})$,

$$\begin{aligned} \langle P_{\mathrm{{DtN}}}^+\phi , \psi \rangle _{\partial {\Omega }}= \langle \phi , P_{\mathrm{{DtN}}}^+\psi \rangle _{\partial {\Omega }}, \end{aligned}$$

where the duality pairing $\langle \phi , \psi \rangle _{\partial {\Omega }}:= \int _{\partial {\Omega }}\phi \, \psi \,\mathrm{d}s$ when $\phi , \psi \in {L^2({\partial {\Omega }})}$; see [19, Equation 2.65, Equation 2.84, Equation A.24]. Therefore, taking the adjoint of (5.1), we find that

$$\begin{aligned} I = \frac{1}{2}I + D'_k- P_{\mathrm{{DtN}}}^+ S_k. \end{aligned}$$

(5.2)

Comparing (5.2) to the definition of $A'_{k,\eta }$ in (1.4), we see that the ideal $\mathrm{i}\eta$ should approximate $P_{\mathrm{{DtN}}}^+$. The idea of choosing $\eta$ as an operator, based on the relations (5.1)-(5.2) (and their analogues for the Neumann problem), essentially first appeared in [55, 56]. The relation (5.1) appeared explicitly in [3, Theorem 2.1], with this paper considering local approximations of the non-local operators $P_{\mathrm{{DtN}}}^+$ and $P_{\mathrm{{NtD}}}^+$, whilst [55, 56] used non-local pseudodifferential-operator approximations.

Regarding 2: In the case when $\partial \Omega$ is a circle (of radius 1), the DtN map is given by

$$\begin{aligned} \frac{\partial u}{\partial r}(1,\theta ) = k\sum _{n=-\infty }^\infty \frac{H_{n}^{(1)'}(k )}{H_{n}^{(1)}(k )} \mathrm{e}^{\mathrm{i}n \theta } \, \left( \frac{1}{2\pi } \int ^{2\pi }_{0}\mathrm{e}^{-\mathrm{i}n \phi }\, u(1,\phi ) \, \mathrm{d}\phi \right) . \end{aligned}$$

The uniform- and double-asymptotic expansions of the Hankel functions (see, e.g., [67, Sections 10.20, 10.41(v)]) imply that

$$\begin{aligned} k \frac{H_{n}^{(1)'}(k )}{H_{n}^{(1)}(k )} \sim {\left\{ \begin{array}{ll} \mathrm{i}k, &{}\text { for }n \text { fixed as } k \rightarrow \infty ,\\ \mathrm{i}k \sqrt{1- \left( \frac{n}{k}\right) ^2}, &{} \,\, n, k \rightarrow \infty \text { with } k-|n| \gg k^{1/3},\\ \mathrm{e}^{2\pi \mathrm{i}/3}\sqrt{\frac{n^2-k^2}{n^{2/3}\zeta }} \frac{ \mathrm{Ai}' \big (\mathrm{e}^{2\pi \mathrm{i}/3}n^{2/3}\zeta \big )}{\mathrm{Ai}\big (\mathrm{e}^{2\pi \mathrm{i}/3}n^{2/3}\zeta \big )},&{} \,\, n, k \rightarrow \infty \text { with } \big ||n|-k\big |\le Mk^{1/3},\\ n\sqrt{1-\left( \frac{k}{n}\right) ^2}, &{} \,\, n, k \rightarrow \infty \text { with } |n|-k\gg k^{1/3}, \end{array}\right. }\nonumber \\ \end{aligned}$$

(5.3)

where $\zeta$ is defined in terms of n and k by [67, Equations 10.20.2 and 10.20.3].³ We see that the approximation $k H_{n}^{(1)'}(k )/H_{n}^{(1)}(k )\sim \mathrm{i}k$ describes the DtN map on the low frequency modes and, in particular, is much better than the approximation $k H_{n}^{(1)'}(k )/H_{n}^{(1)}(k )\sim -\mathrm{i}k$. The asymptotics (5.3), however, show that neither the approximations $\mathrm{i}k$ or $-\mathrm{i}k$ are particularly good on the higher frequency modes. An almost-identical analysis is valid for the sphere, and more generally for a smooth convex curved obstacle, since the symbol of the DtN map for such domains is described by the asymptotics (5.3); see [37, Section 9, last formula on p. 58].