A high frequency boundary element method for scattering by a class of multiple obstacles

TL;DR

This paper introduces a high-frequency boundary element method combining hybrid numerical asymptotics for convex polygons with standard methods for other obstacles, significantly reducing degrees of freedom needed at high frequencies.

Contribution

The paper develops a novel boundary element method that achieves logarithmic growth in degrees of freedom with frequency for scattering problems involving multiple obstacles, especially with large convex polygons.

Findings

Degrees of freedom grow logarithmically with frequency for the proposed method.

The method is most effective for large convex polygons and small obstacles with combined perimeter comparable to the wavelength.

Significant reduction in computational complexity compared to standard polynomial schemes.

Abstract

We propose a boundary element method for problems of time-harmonic acoustic scattering by multiple obstacles in two dimensions, at least one of which is a convex polygon. By combining a Hybrid Numerical Asymptotic (HNA) approximation space on the convex polygon with standard polynomial-based approximation spaces on each of the other obstacles, we show that the number of degrees of freedom required in the HNA space to maintain a given accuracy needs to grow only logarithmically with respect to the frequency, as opposed to the (at least) linear growth required by standard polynomial-based schemes. This method is thus most effective when the convex polygon is many wavelengths in diameter and the small obstacles have a combined perimeter comparable to the problem wavelength.

Equations

$$u^i(\mathbf{x}):={\mathrm{e}}^{\mathrm{i}k\mathbf{x}\cdot \mathbf{d}},\mathbf{x}\in {\mathbb{R}}^2,$$

$$\Delta u+k^{2}u=0$$

$$u=0$$

$$\left(\frac{\partial }{\partial r}-\mathrm{i}k\right)u^s(\mathbf{x})=o(r^{-1/2}),\text{as }r:=|\mathbf{x}|\to \infty .$$

$$S_k\phi (\mathbf{x}):={\int }_{\partial D}{\Phi }_k(\text{x},\text{y})\phi (\text{y})ds(\text{y}),\mathbf{x}\in D,$$

$$u=u^i-S_k\frac{\partial u}{\partial \mathbf{n}}\text{in }D.$$

$${\mathcal{A}}_{k,\eta }:=\frac{1}{2}\mathcal{I}+{\mathcal{D}}_k^{\prime }-\mathrm{i}\eta {\mathcal{S}}_k,$$

$${\mathcal{S}}_k\phi (\mathbf{x}):={\int }_{\partial D}{\Phi }_k(\mathbf{x},\mathbf{y})\phi (\mathbf{y})ds(\mathbf{y}),\mathbf{x}\in \partial D,\phi \in L^2(\partial D),$$

$${\mathcal{D}}_k^{\prime }\phi (\mathbf{x}):={\int }_{\partial D}\frac{\partial {\Phi }_k(\mathbf{x},\mathbf{y})}{\partial \mathbf{n}(\mathbf{x})}\phi (\mathbf{y})ds(\mathbf{y}),\mathbf{x}\in \partial D,\phi \in L^2(\partial D),$$

$${\mathcal{A}}_{k,\eta }\frac{\partial u}{\partial \mathbf{n}}=f_{k,\eta },\text{on }\partial D,$$

$$f_{k,\eta }=\left(\frac{\partial }{\partial \mathbf{n}}-\mathrm{i}\eta \right)u^i.$$

$${\left({\mathcal{A}}_{k,\eta }\frac{\partial u}{\partial \mathbf{n}},w\right)}_{L^2(\Gamma \cup \gamma )}={\left(f_{k,\eta },w\right)}_{L^2(\Gamma \cup \gamma )},\text{for all }w\in L^2(\Gamma \cup \gamma ),$$

$$\mathbf{Z}(\mathbf{x}):=(x_1\chi (\mathbf{x}),x_2),\mathbf{x}=(x_1,x_2)\in {\mathbb{R}}^2.$$

$$G_j(\mathbf{x},\mathbf{y}):={\Phi }_k(\mathbf{x},\mathbf{y})-{\Phi }_k({\tilde{\mathbf{x}}}^j,\mathbf{y}),\mathbf{x}\ne \mathbf{y},$$

$$\frac{\partial G_j(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_j(\mathbf{y})}=2\frac{\partial {\Phi }_k(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_j(\mathbf{y})}\text{and}{G}_j(\mathbf{x},\mathbf{y})=0,\text{for}\mathbf{y}\in {\Gamma }_j^{\infty }.$$

$$u^s(\mathbf{x})=$$

$$+{\int }_{\gamma \cap U_j}\left[\frac{\partial G_j(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_{\gamma }(\mathbf{y})}{u}^s(\mathbf{y})-G_j(\mathbf{x},\mathbf{y})\frac{\partial u^s(\mathbf{y})}{\partial {\mathbf{n}}_{\gamma }}\right]ds(\mathbf{y})$$

$$-{\int }_{\partial B_R\cap U_j}\left[\frac{\partial G_j(\mathbf{x},\mathbf{y})}{\partial r}{u}^s(\mathbf{y})-G_j(\mathbf{x},\mathbf{y})\frac{\partial u^s(\mathbf{y})}{\partial r}\right]ds(\mathbf{y}),$$

$$\displaystyle u^{s}(\mathbf{x})=2\int_{\Gamma^{\infty}_{j}\setminus{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}}\frac{\partial\Phi_{k}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{j}(\mathbf{y})}u^{s}(\mathbf{y})\operatorname{d}\!{s(\mathbf{y})}+\int_{\gamma\cap U_{j}}\left[\frac{\partial G_{j}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{\gamma}(\mathbf{y})}u^{s}(\mathbf{y})-G_{j}(\mathbf{x},\mathbf{y})\frac{\partial u^{s}(\mathbf{y})}{\partial\mathbf{n}_{\gamma}}\right]\operatorname{d}\!{s(\mathbf{y})},$$

$$\displaystyle\Bigg{(}\int_{\gamma\cap U_{j}}-$$

$$\displaystyle=\int_{U_{j}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}}\left[\Delta G_{j}(\mathbf{x},\mathbf{y})u^{i}(\mathbf{y})-G_{j}(\mathbf{x},\mathbf{y})\Delta u^{i}(\mathbf{y})\right]\operatorname{d}\!V(\mathbf{y})=0,$$

$${\int }_{\gamma \cap U_j}\left[\frac{\partial G_j(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_{\gamma }(\mathbf{y})}{u}^s(\mathbf{y})-G_j(\mathbf{x},\mathbf{y})\frac{\partial u^s(\mathbf{y})}{\partial {\mathbf{n}}_{\gamma }}\right]ds(\mathbf{y})$$

$$=$$

$$\stackrel{\text{eqref}G{2}_{i}ngamma}{=}$$

$$u^s(\mathbf{x})=$$

$$\displaystyle-2\int_{\Gamma_{j}^{\infty}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}}\frac{\partial\Phi_{k}(\mathbf{x},\mathbf{y})}{\partial\mathbf{n}_{j}(\mathbf{y})}u^{i}(\mathbf{y})\operatorname{d}\!{s(\mathbf{y})},\qquad\mathbf{x}\in D\cap U_{j}.$$

$$u^i(\mathbf{x})=2{\int }_{{\Gamma }_j^{\infty }}\frac{\partial {\Phi }_k(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_j(\mathbf{y})}{u}^i(\mathbf{y})ds(\mathbf{y}),\mathbf{x}\in U_j.$$

$$\frac{\partial u(\mathbf{x})}{\partial \mathbf{n}}=$$

$$-2{\int }_{\gamma \cap U_j}\frac{\partial {\Phi }_k(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_j(\mathbf{x})}\frac{\partial u(\mathbf{y})}{\partial {\mathbf{n}}_{\gamma }}ds(\mathbf{y}),\mathbf{x}\in {\Gamma }_j,{\mathbf{n}}_j\cdot \mathbf{d}\ge 0.$$

$$u^i(\mathbf{x})=-2{\int }_{{\Gamma }_j^{\infty }}\frac{\partial {\Phi }_k(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_j(\mathbf{y})}{u}^i(\mathbf{y})ds(\mathbf{y}),\mathbf{x}\in {\mathbb{R}}^2\setminus {\overline{U}}_j.$$

$$u^r(\mathbf{x})=2{\int }_{{\Gamma }_j^{\infty }}\frac{\partial {\Phi }_k({\tilde{\mathbf{x}}}^j,\mathbf{y})}{\partial {\mathbf{n}}_j(\mathbf{y})}{u}^i(\mathbf{y})ds(\mathbf{y})=-2{\int }_{{\Gamma }_j^{\infty }}\frac{\partial {\Phi }_k(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_j(\mathbf{y})}{u}^i(\mathbf{y})ds(\mathbf{y}).$$

$$u^i(\mathbf{x})=u^i(\mathbf{x})+u^r(\mathbf{x})+2{\int }_{{\Gamma }_j^{\infty }}\frac{\partial {\Phi }_k(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_j(\mathbf{y})}{u}^i(\mathbf{y})ds(\mathbf{y}).$$

$$\frac{\partial u(\mathbf{x})}{\partial \mathbf{n}}=$$

$$-2{\int }_{\gamma \cap U_j}\frac{\partial {\Phi }_k(\mathbf{x},\mathbf{y})}{\partial {\mathbf{n}}_j(\mathbf{x})}\frac{\partial u(\mathbf{y})}{\partial {\mathbf{n}}_{\gamma }}ds(\mathbf{y}),\mathbf{x}\in {\Gamma }_j,{\mathbf{n}}_j\cdot \mathbf{d}<0,$$

$$\Psi(\mathbf{x}):=\left\{\begin{array}[]{cl}2{\partial u^{i}(\mathbf{x})}/{\partial\mathbf{n}},&\mathbf{x}\in\Gamma_{j}\subset\Gamma:\mathbf{n}_{j}(\mathbf{x})\cdot\mathbf{d}<0,\\ 0,&\mathbf{x}\in\Gamma_{j}\subset\Gamma:\mathbf{n}_{j}(\mathbf{x})\cdot\mathbf{d}\geq 0.\end{array}\right.$$

Full text

Helmholtz, high frequency, multiple scattering, integral equations, BEM, hp discretisation, HNA method.

1 Introduction

Standard finite or boundary element methods for wave scattering problems, with piecewise polynomial approximation spaces, typically require at least a fixed number of degrees of freedom per wavelength to maintain accuracy as the frequency of the incident wave increases. This dependence can lead to a requirement for an excessively large number of degrees of freedom at high frequencies.

For certain geometries, the Hybrid Numerical Asymptotic (HNA) approach (see, e.g., Chandler-Wilde et al. (2012b) and the references therein) overcomes this restriction by absorbing the high frequency asymptotic behaviour into the approximation space. This is implemented via a Boundary Element Method (BEM), which is particularly effective as the high frequency behaviour need only be captured on the surface of the obstacle. When constructing an HNA method, a key ingredient is an understanding of the high frequency asymptotics of the underlying physical problem. Such methods are well-studied for strictly convex smooth scatterers, see for example Bruno et al. (2004); Domínguez et al. (2007); Ecevit & Özen (2017); Ecevit (2018); Ecevit & Eruslu (2019), and the many references therein. This paper is instead concerned with HNA for polygonal scatterers, specifically extending this method to multiple obstacles. We discuss the relevant literature in detail now.

In contrast with standard BEM, in which the number of degrees of freedom (DOFs) required to accurately represent the solution depends linearly on the frequency, the number of DOFs needed to achieve a given accuracy (for scattering by a convex polygon in two-dimensions) was shown to depend only logarithmically on the frequency for the $h$-BEM version of HNA in Chandler-Wilde & Langdon (2007), this improved to the $hp$-BEM version in Hewett et al. (2013). These ideas were extended, in Chandler-Wilde et al. (2015), to a certain class of non-convex polygons, with the high frequency asymptotics arising from re-reflections and partial illumination (shadowing) being fully captured by a careful choice of approximation space. Similar ideas have been applied to penetrable obstacles in Groth et al. (2015, 2018) and to two- and three-dimensional screens in Hewett et al. (2015) and Hargreaves et al. (2015) respectively. All of these methods are, broadly speaking, for single obstacles and for plane wave incidence (although an extension to other incident fields is discussed in Remark 3.14). In this paper we extend the HNA method to a class of more general multiple scattering configurations. Although here we focus on the case where at least one of the obstacles is a convex polygon, the ideas we present may be applied to the same problem where this obstacle is (for example) strictly convex and smooth. The key ingredient is that this obstacle be amenable, for the corresponding single scattering problem, to solution by an HNA-BEM.

Problems of high frequency scattering by one large relatively simple obstacle and one (or many) small obstacle(s) are potentially of practical interest. An approach used in Lenoir et al. (2017) and Bonnet et al. (2018) for such problems is to appeal to high frequency asymptotics on the large obstacle, via Geometrical/Physical Optics approximation, and approximate the solution on (or in some neighbourhood of) the small components using a standard BEM/FEM. This approximation works well at sufficiently high frequencies, but ignores diffracted waves emanating from the large obstacle, and so is not controllably accurate across all frequencies. Moreover, a Geometrical Optics approach will include a ray-tracing algorithm, which typically requires that the multiple scattering problem is solved iteratively, see Ecevit & Reitich (2009); Anand et al. (2010); Geuzaine et al. (2005). This involves reformulation as a Neumann series consisting entirely of operators on a single scatterer. More generally, iterative approaches are common in multiple scattering problems and work well for certain configurations. However, the iterative approach cannot be applied to all such problems: the Neumann series will diverge for cases where the separation of the obstacles is too small. A method which accelerates the rate of convergence in the Neumann series is presented in Boubendir et al. (2017).

In this paper, we present a method which is particularly effective for high-frequency time-harmonic scattering by one large obstacle and one (or many) small obstacle(s). Specifically, since we use an oscillatory basis on the large obstacle, this can be many wavelengths long. On the small obstacle(s), we use a piecewise-polynomial basis, so the method is most effective when their combined perimeter is comparable to one wavelength. In contrast to the other methods currently available for similar problems, the method we present in this paper is controllably accurate and does not need to be solved iteratively, whilst the only constraint on the separation of the obstacles is that they must be $O(\lambda)$ apart, where $\lambda$ denotes wavelength, hence the obstacles may be very close together at high frequencies.

1.1 Outline of the paper

In §2 we describe in detail the class of multiple scattering problem we are aiming to solve.

In §3.1 we extend the representation on which the HNA-BEM for a single convex polygon is based (Chandler-Wilde & Langdon, 2007; Hewett et al., 2011), to account for the contribution to the solution from neighbouring scatterers. In particular, we derive a representation (3.1), which decomposes the Neumann trace on the large obstacle into the sum of a physical optics term, the diffracted waves, and the contribution arising from interaction with the small obstacles. In Theorem 3.6 we bound the Helmholtz solution in the complement of the scatterers with bounding constant algebraic in the wavenumber. We describe the singular behaviour of the envelope of the diffracted waves in Theorem 3.10.

In §4.1 we construct an approximation space enriched with oscillatory basis functions, designed to represent the solution on the large obstacle with a number of degrees of freedom that does not need to increase significantly as frequency grows. We also describe some potential advantages of using an approximation space based on a single mesh, as opposed to the overlapping meshes of Hewett et al. (2015); Chandler-Wilde & Langdon (2007); Hewett et al. (2013). In addition, in §4.2 we define a (standard) piecewise-polynomial space on the small scatterer. Conditions for exponential convergence of the best approximation in terms of polynomial degree on the large and small scatterers are given by Theorem 4.3 and Proposition 4.7 respectively.

In §5 a Galerkin method using this approximation space is outlined, alongside related error estimates of the total field and far-field pattern. Numerical results for an implementation of this method are presented in §6. While Theorem 5.3 ensures exponential convergence when the big scatterer is a convex polygon and the small scatterers are analytic, §6 demonstrates that exponential convergence is still possible when the small scatterers have corners.

Finally, in the appendix, we introduce an alternative boundary integral equation formulation, which is provably coercive under certain geometric constraints. This gives us explicit quasi-optimality estimates, which when combined with results in earlier sections could be used to give explicit error estimates for a certain class of multiple scattering configurations.

2 Problem statement

We consider the two-dimensional problem of time-harmonic acoustic scattering by ${\mathcal{N}_{\gamma}}+1$ sound-soft scatterers, at least one of which is an ${\mathcal{N}_{\Gamma}}$-sided convex polygon. In addition to this convex polygon, we assume that the other ${\mathcal{N}_{\gamma}}$ obstacles are pairwise disjoint with Lipschitz piecewise-$C^{1}$ boundary. Denote the interior of the convex polygon by {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\subset\mathbb{R}^{2} and its boundary by \Gamma:=\partial{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}. We denote by $\Gamma_{j}$ the $j$th side of $\Gamma$, for $j=1,\ldots,{\mathcal{N}_{\Gamma}}$. The bounded open set {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}:=\bigcup_{i=1}^{{\mathcal{N}_{\gamma}}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{i}}\subset\mathbb{R}^{2}\setminus\overline{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}}} represents the collection of the ${\mathcal{N}_{\gamma}}$ other obstacles, which are denoted {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{i}}, for $i=1,\ldots,{\mathcal{N}_{\gamma}}$. We denote the combined Lipschitz boundary of these by \gamma:=\partial{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}. The unbounded exterior domain is denoted D:=\mathbb{R}^{2}\setminus(\overline{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}}\cup\overline{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}}), with boundary $\partial D=\Gamma\cup\gamma$. The normal derivative operator (or Neumann trace) is defined as $\partial/\partial\mathbf{n}:=\mathbf{n}\cdot\nabla$, in which $\mathbf{n}=(n_{1},n_{2})$ denotes the unit normal directed into $D$; we denote also $\mathbf{n}_{j}:=\mathbf{n}|_{\Gamma_{j}}$ and $\mathbf{n}_{\gamma}:=\mathbf{n}|_{\gamma}$. We assume that the distance between {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} is positive, so that ${\partial D}$ is a Lipschitz boundary. A simple example of a geometric configuration that fits inside of this framework is depicted in Figure 1. We note that throughout the paper, it is the quantities $\Gamma$ and $\gamma$ which are used most frequently.

We aim to solve the following boundary value problem (BVP): given the incident plane wave

[TABLE]

where $k:=2\pi/\lambda>0$ denotes the wavenumber (for wavelength $\lambda$) and $\mathbf{d}\in\mathbb{R}^{2}$ is a unit direction vector, determine the total field $u\in C^{2}(D)\cap C(\bar{D})$ such that

[TABLE]

and $u^{s}:=u-u^{i}$ satisfies the Sommerfeld radiation condition Colton & Kress (2013, (3.62))

[TABLE]

Problems for a broader class of incident field $u^{i}$ are discussed briefly in Remark 3.14.

The BVP (2)–(4) can be reformulated as a boundary integral equation (BIE). We denote the single layer potential $S_{k}:L^{2}(\partial D)\rightarrow C^{2}(D)$ by

[TABLE]

where $\Phi_{k}(\mathbf{x},\mathbf{y}):=(\mathrm{i}/4)H^{(1)}_{0}(k|\mathbf{x}-\mathbf{y}|)$ is the fundamental solution of (2), in which $H^{(1)}_{0}$ denotes the Hankel function of the first kind and order zero. If $u$ satisfies the BVP (2)–(4), then ${\partial u}/{\partial\mathbf{n}}\in L^{2}(\partial D)$ and the following Green’s representation holds (see, e.g., Chandler-Wilde et al. (2012b, Theorem 2.43))

[TABLE]

Definition 2.1** (Combined potential operator).**

The standard combined potential operator ${\mathcal{A}}_{k,\eta}:L^{2}(\partial D)\rightarrow L^{2}(\partial D)$ (see, e.g., Colton & Kress (2013); Chandler-Wilde et al. (2012b)) is defined by

[TABLE]

where ${\mathcal{I}}$ is the identity operator, $\eta\in\mathbb{R}\setminus\{0\}$ is a coupling parameter,

[TABLE]

denotes the single layer operator and

[TABLE]

denotes the adjoint of the double-layer operator.

From (6), the BVP (2)–(3) can be reformulated as a BIE (Chandler-Wilde et al., 2012b, (2.69), (2.114))

[TABLE]

where the right-hand side data $f_{k,\eta}\in L^{2}(\partial D)$ is

[TABLE]

It follows from Chandler-Wilde et al. (2012b, Theorem 2.27) that ${\mathcal{A}}_{k,\eta}$ is invertible. A standard variational form of (7) is

[TABLE]

which can be approximated by a piecewise-polynomial Galerkin BEM. Our approach differs from this in that: (i) we decompose the unknown ${\partial u}/{\partial\mathbf{n}}$ into a known physical optics term, a diffracted term, and a term which expresses the leading order behaviour on $\Gamma$ in terms of the solution on $\gamma$ (see §3.1); (ii) we approximate the diffracted term on $\Gamma$ using an oscillatory basis (see §4.1). The use of this basis is justified by the representation and regularity results in §3.3. This leads to a new variational formulation that is equivalent to (9). This is shown in (57)–(58) with the resulting Galerkin scheme given in equations (59)–(60).

2.1 Geometric assumptions

In related literature, there appears to be no single consistent definition of the term polygon, so we shall clarify a definition that is appropriate for this paper.

Definition 2.2** (Polygon).**

We say $\Upsilon{\subset\mathbb{R}^{2}}$ is a polygon if it is a bounded Lipschitz open set with a boundary consisting of a finite number of straight line segments.

We note that Definition 2.2 permits multiple disconnected shapes, whereas other conventions in related literature do not. As we impose that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} is convex, it cannot consist of disconnected components. On the other hand, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} may consist of disconnected components. Many results that follow hold for a subclass of polygons, which we define now (as in, e.g., Spence (2014, Definition 1.1)).

Definition 2.3** (Non-trapping polygon).**

We say that a polygon $\Upsilon$ (in the sense of Definition 2.2) is non-trapping if:

1. (i)

No three vertices of $\Upsilon$ are co-linear; 2. (ii)

For a ball $B_{R}$ with radius $R>0$ sufficiently large that $\Upsilon\subset B_{R}$, there exists a $T(R)<\infty$ such that all billiard trajectories that start inside of $B_{R}\setminus\Upsilon$ at time $T=0$ and miss the vertices of $\Upsilon$ will leave $B_{R}$ by time $T(R)$.

Previous analyses of HNA methods (e.g., Hewett et al. (2013); Chandler-Wilde et al. (2015)) have instead relied upon convergence and regularity estimates for scattering obstacles which are convex or star-shaped (introduced formally in Definition A.2), a property not enjoyed by multiple scattering configurations. However, for configurations which satisfy the conditions of Definition 2.3, bounds on the Dirichlet-to-Neumann (DtN) maps are known (Baskin et al., 2016), which will provide an alternative route to bounding the solution to (2)–(4) in §3.2.

In addition to the theory of Baskin et al. (2016) for non-trapping polygons, we shall consider a certain class of trapping configurations, for which bounds on DtN maps were recently derived in Chandler-Wilde et al. (2020), building on the earlier work of Galkowski et al. (2019); Baskin et al. (2016). These estimates will form a key component of our numerical analysis, in particular enabling us to bound the solution to (2)–(4) in §3.2, and obtain best approximation on $\gamma$ in Proposition 4.7. A formal definition of these so-called $(R_{0},R_{1})$ configurations will follow, but these may be loosely interpreted as configurations $\Upsilon$ which are star-shaped outside of some ball. There is a second ball inside of the first, whose radius is sufficiently small, and inside of which some trapping may occur.

Definition 2.4** ($(R_{0},R_{1})$ configuration).**

For $0<R_{0}<R_{1}$ we say that a Lipschitz $\Upsilon$ is an $(R_{0},R_{1})$ configuration if there exists a $\chi\in C^{3}[0,\infty)$ which satisfies

1. (i)

$\chi(|\mathbf{x}|)=0$* for $0\leq|\mathbf{x}|\leq R_{0}$, $\chi(|\mathbf{x}|)=1$ for $|\mathbf{x}|\geq R_{1},0<\chi(|\mathbf{x}|)<1,$ for $R_{0}\leq|\mathbf{x}|\leq R_{1}$,* 2. (ii)

$0\leq\chi^{\prime}(|\mathbf{x}|)\leq 4$, for $|\mathbf{x}|>0$,

such that $\mathbf{Z}(\mathbf{x})\cdot\mathbf{n}(\mathbf{x})\geq 0$ for all $\mathbf{x}\in\partial\Upsilon$ for which the normal $\mathbf{n}(\mathbf{x})$ is defined, where

[TABLE]

Naturally, one can rotate the coordinate system if required to ensure the above conditions hold. For further explanation and examples of $(R_{0},R_{1})$ configurations, we refer to Chandler-Wilde et al. (2020, §1.2.1).

3 Representation and regularity of solution on $\Gamma$

The structure of this section is as follows: In §3.1 we extend the single scattering HNA ansatz of Chandler-Wilde & Langdon (2007, (3.5)) to a multiple scattering configuration, introducing a new operator which accounts for the other obstacle(s). In §3.2 we bound the solution of the multiple scattering problem in the domain $D$, a necessary component of the best approximation estimates that follow. In §3.3 we show that the envelopes of the diffracted waves, which the HNA space is designed to approximate, behave similarly to the single scattering problem (under very reasonable assumptions). This means that the HNA space of Hewett et al. (2013) may be used on the convex polygon in the multiple scattering approximation without any modification (though with a different leading order term).

3.1 The representation formula for the Neumann trace on $\Gamma$

As in Chandler-Wilde & Langdon (2007, §3), we will extend a single side $\Gamma_{j}$ of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}, and solve the resulting half-plane problem, to obtain an explicit representation for ${\partial u}/{\partial\mathbf{n}}$ on $\Gamma_{j}$ in terms of known oscillatory functions on $\Gamma_{j}$ and (in contrast to Chandler-Wilde & Langdon (2007, §3)) ${\partial u}/{\partial\mathbf{n}}$ on $\gamma$. This representation will form the ansatz used for the discretisation. Throughout this section, when $u$ or $u^{s}$ is restricted to $\Gamma\cup\gamma$, it is assumed that the exterior trace has been taken. Considering a single side $\Gamma_{j}$ of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}, $1\leq j\leq{\mathcal{N}_{\Gamma}}$, define $\Gamma_{j}^{+}$ and $\Gamma_{j}^{-}$ as the infinite extensions of $\Gamma_{j}$, each as a straight half line in the clockwise and anti-clockwise direction (about the interior {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}) respectively (see Figure 2). Denote by $U_{j}$ the (open) half-plane with boundary $\Gamma_{j}^{\infty}:=\Gamma_{j}^{-}\cup\Gamma_{j}\cup\Gamma_{j}^{+}$, chosen such that $U_{j}$ does not contain {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}. We informally call $U_{j}$ the upper half-plane relative to $\Gamma_{j}$. On $\Gamma_{j}^{\infty}$, the unit normal $\mathbf{n}_{j}$ points into $U_{j}$. Define the half-plane Dirichlet Green’s function

[TABLE]

where ${{\widetilde{\mathbf{x}}}^{j}}$ is the reflection of $\mathbf{x}$ across $\Gamma_{j}$. Formally, $\mathbf{x}={{\widetilde{\mathbf{x}}}^{j}}$ when $\mathbf{x}\in\Gamma_{j}$, otherwise ${{\widetilde{\mathbf{x}}}^{j}}\neq\mathbf{x}$ satisfies $\text{dist}(\{\mathbf{x}\},\Gamma_{j}^{\infty})=\text{dist}(\{{{\widetilde{\mathbf{x}}}^{j}}\},\Gamma_{j}^{\infty})=\tfrac{1}{2}|\mathbf{x}-{{\widetilde{\mathbf{x}}}^{j}}|$. It follows that

[TABLE]

We let $B_{R}$ be an open ball of radius $R$ centred at the origin, with $R$ chosen sufficiently large that U_{j}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}\subset B_{R}, i.e. all the scatterers in the relative upper half-plane lie inside the ball.

Green’s second identity can now be applied to $G_{j}(\mathbf{x},\cdot)$ and $u^{s}{=u-u^{i}}$ for ${\mathbf{x}~{}\in}~{}D~{}\cap~{}U_{j}~{}\cap~{}B_{R}$, together with (10), to obtain

[TABLE]

where ${\partial}/{\partial\mathbf{n}_{j}}=\mathbf{n}_{j}\cdot\nabla$ and ${\partial}/{\partial\mathbf{n}_{\gamma}}=\mathbf{n}_{\gamma}\cdot\nabla$, $\mathbf{n}_{j}$ and $\mathbf{n}_{\gamma}$ are the unit normal vector fields pointing into $D\cap U_{j}\cap B_{R}$ from \Gamma^{\infty}_{j}\cap B_{R}\setminus{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} and from $\gamma\cap U_{j}$, respectively, and ${\partial}/{\partial r}=\frac{\mathbf{y}}{|\mathbf{y}|}\cdot\nabla$ denotes the normal derivative on $\partial B_{R}\cap U_{j}$ pointing out of $D\cap U_{j}\cap B_{R}$. As $R\rightarrow\infty$, the third integral vanishes by the same reasoning as in, e.g., Colton & Kress (2013, Theorem 2.4). The representation (11) then becomes

[TABLE]

for \mathbf{x}\in U_{j}\setminus{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}.

We now apply Green’s second identity to $u^{i}$ and $G_{j}(\mathbf{x},\mathbf{y})$ in U_{j}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} and obtain, for $\mathbf{x}\in D\cap U_{j}$,

[TABLE]

as $u^{i}$ and $\Phi_{k}(\mathbf{x},\cdot)$ satisfy the Helmholtz equation (2) in {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} for $\mathbf{x}\in D\cap U_{j}$. The sign of the boundary integral differs on the two parts of \partial(U_{j}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega})=(\gamma\cap U_{j})\cup(\Gamma_{j}^{\infty}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}) because the normal derivative ${\partial}/{\partial\mathbf{n}}$ involves the outward-pointing normal vector $\mathbf{n}_{\gamma}$ on $\gamma\cap U_{j}$ and the inward-pointing normal $\mathbf{n}_{j}$ on \Gamma_{j}^{\infty}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}, as depicted in Figure 2.

We then use $u^{s}=u-u^{i}$ to expand the last term in (12): for $\mathbf{x}\in D\cap U_{j}$

[TABLE]

Substituting this expression in (12) and using again (10), we obtain a representation for $u^{s}$:

[TABLE]

The final term will be non-zero only if \Gamma_{j}^{\infty}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}\neq\emptyset, namely, in case one of the components of $\gamma$ is $\Gamma_{j}^{\infty}$ (see e.g. the component {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{3}} in Figure 2).

This integral representation must be combined with one for $u^{i}$ to construct a useful representation for ${\partial u}/{\partial\mathbf{n}}$ on $\Gamma$. The half-plane representation of Chandler-Wilde (1997, §3) can be applied to upward propagating plane waves. We consider first the case $\mathbf{n}_{j}\cdot\mathbf{d}\geq 0$, which means that $\Gamma_{j}$ is in shadow, from Chandler-Wilde & Langdon (2007, (3.3)):

[TABLE]

Adding this to (14) and taking the Neumann trace on $\Gamma_{j}$, we obtain a representation for the solution

[TABLE]

For a downward-propagating wave $\mathbf{n}_{j}\cdot\mathbf{d}<0$, i.e. when $\Gamma_{j}$ is illuminated by $u^{i}$, we can apply the same result to the lower half-plane $\mathbb{R}^{2}\setminus\overline{U}_{j}$ (where the direction of the normal is reversed)

[TABLE]

Now define $u^{r}(\mathbf{x}):=-u^{i}({\widetilde{\mathbf{x}}}^{j})$ for $\mathbf{x}\in U_{j}$. Intuitively, $u^{r}$ may be considered the reflection of $u^{i}$ by a sound-soft line at $\Gamma_{j}^{\infty}$. It follows that ${\partial u^{r}}/{\partial\mathbf{n}_{j}}={\partial u^{i}}/{\partial\mathbf{n}_{j}}$ on $\Gamma_{j}^{\infty}$ and, for $\mathbf{x}\in U_{j}$,

[TABLE]

Rearranging this and adding $u^{i}$ gives

[TABLE]

Summing with (14) and taking the Neumann trace gives the representation for ${\partial u}/{\partial\mathbf{n}}$ on $\Gamma_{j}$:

[TABLE]

where we used again (10) and ${\partial u^{r}}/{\partial\mathbf{n}_{j}}={\partial u^{i}}/{\partial\mathbf{n}_{j}}$ on $\Gamma_{j}$.

The representation (15)–(16) may be viewed as a correction to the Physical Optics approximation for a single scatterer, which is defined as

[TABLE]

Specifically, this correction can be split into two parts. The first integral of (15) and (16) represents the waves diffracted by the corners of $\Gamma$ (diffraction is ignored by the Physical Optics approximation ), whilst the second integral represents the correction to the waves reflected by the sides of $\Gamma$, as a result of the presence of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}. Unless the distance between the scatterers is sufficiently large, it is reasonable to expect the second correcting term to be not negligible; see Gibbs (2017, Lemma 4.2) for a precise quantification of this fact.

We now write more explicitly the integral representation (15)–(16) in terms of the parametrisations of the segments $\Gamma_{j}$ and of their extensions $\Gamma_{j}^{\infty}$. From the standard properties of Bessel functions (see, e.g., DLMF (2019, §10)), we have that for $\mathbf{x}\in\Gamma_{j}$, \mathbf{y}\in\Gamma^{\pm}_{j}\setminus{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega},

[TABLE]

see Chandler-Wilde & Langdon (2007, (3.6)). To make use of this identity, we parametrise $\Gamma$ by

[TABLE]

where $L_{j}$ is the length of the $j$th side, $\mathbf{P}_{j}$ is the $j$th corner of $\Gamma$, and ${\widetilde{L}}_{j}:=\sum_{\ell=1}^{j}{L}_{\ell}$ is the arc length up to the $(j+1)$th corner, with $\mathbf{P}_{{{\mathcal{N}_{\Gamma}}}+1}:=\mathbf{P}_{1}$. We will also denote by $L_{\Gamma}:={\widetilde{L}}_{{\mathcal{N}_{\Gamma}}}$ the total length of $\Gamma$. Similarly we parametrise $\Gamma_{j}^{-}\cup\Gamma_{j}\cup\Gamma_{j}^{+}$ by

[TABLE]

We use (15)–(16) to represent the solution on a single side $\Gamma_{j}$, extending the ansatz of Chandler-Wilde & Langdon (2007); Hewett et al. (2013) to multiple scattering problems

[TABLE]

we shall now discuss each term in the ansatz separately. Here $\Psi$ is the Physical Optics approximation (17), with the envelopes of the diffracted waves on each side defined by

[TABLE]

where $Z_{j}^{+}:=\{t\in\mathbb{R}:\mathbf{y}_{j}({\widetilde{L}}_{j-1}-t)\in\gamma\}$ and $Z_{j}^{-}:=\{t\in\mathbb{R}:\mathbf{y}_{j}({\widetilde{L}}_{j}+t)\in\gamma\}$ are used to exclude from the integral the points inside {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} (as is the case for {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{3}} of Figure 2), to remain consistent with (15)–(16). The interaction operator ${\mathcal{G}}_{{\gamma\shortrightarrow\Gamma}_{j}}:L^{2}(\gamma)\to L^{2}(\Gamma_{j})$ used in (3.1) is based on the final term of (15)–(16), and is defined by

[TABLE]

for $\varphi\in L^{2}(\gamma)$. We extend this definition to ${\mathcal{G}}_{\gamma\shortrightarrow\Gamma}:L^{2}(\gamma)\to L^{2}(\Gamma)$ as

[TABLE]

Remark 3.1**.**

The ansatz (3.1) is an extension of Chandler-Wilde & Langdon (2007, (3.9)) and Hewett et al. (2013, (3.2)), with an additional term which relates the solution on $\Gamma$ to the solution on $\gamma$. It is important to note that this additional term is not the only term influenced by the presence of $\gamma$ and that one cannot solve for $v_{\pm}$ on a single scatterer and then add the ${\mathcal{G}}_{\gamma\shortrightarrow\Gamma}[{\partial u}/{\partial\mathbf{n}|_{\gamma}}]$ term. The reason for this is clear from (20)–(21): even if $Z_{j}^{\pm}$ were of measure zero, so that the equations for (20)–(21) were identical to the case of a single scatterer, the integral contains $u$, which depends on the configuration $\partial D$. Intuitively this makes sense, diffracted waves emanating from the corners of $\Gamma$ will also be influenced by the presence of additional scatterers.

Many of the bounds which follow are explicit only in $k$ or the parameters which determine meshwidth or polynomial degree of an approximation space. Henceforth we will use $A\lesssim B$ to mean $A\leq CB$, where $C$ is a constant that depends only on the geometry of $\Upsilon$. To gauge the size of the contribution to the reflected waves on $\Gamma$ arising from the presence of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}, we require the following bound on the operator ${\mathcal{G}}_{\gamma\shortrightarrow\Gamma}$.

Lemma 3.2**.**

For $\partial D=\Gamma\cup\gamma$ with $\Gamma$ and $\gamma$ disjoint, we have the following bound on the interaction operator ${\mathcal{G}}_{\gamma\shortrightarrow\Gamma}$ defined in (23), given $k_{0}>0$:

[TABLE]

where

[TABLE]

where $L_{\Gamma}$ and $L_{\gamma}$ denote the perimeters of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} respectively.

Proof 3.3**.**

For $0\neq\varphi\in L^{2}(\gamma)$, using the Cauchy–Schwarz inequality, we can write

[TABLE]

The result follows from ${H^{(1)}_{0}}^{\prime}(z)=-H^{(1)}_{1}(z)$ and Chandler-Wilde et al. (2009, (1.23)), which states that $|H^{(1)}_{1}(z)|$ $\leq\sqrt{2/(\pi z)}+{2/(\pi z)}$ for $z>0$.

As intuition would suggest, Lemma 3.2 confirms that the norm of the interaction operator (23) decreases as the obstacles move further apart, i.e., as the interaction between them decreases.

3.2 Estimates of the $L^{\infty}$ norm of the Helmholtz solution in $D$

A value that will feature in many of the estimates for this method is

[TABLE]

The dependence of $u_{\max}(k)$ on the wavenumber $k$ is of key importance, as $u_{\max}(k)$ appears as a multiplicative constant in the $hp$ best approximation result derived in §5, alongside a term which decreases exponentially with $p$. To show exponential convergence of the method, we therefore require that $u_{\max}(k)$ grows at most algebraically with $k$. To explore this dependence, we will make use of the current best available bounds on the Dirichlet-to-Neumann map (see, e.g., Chandler-Wilde et al. (2012b, §2.7)) for multiple obstacle configurations. Recently in Chandler-Wilde et al. (2020) such bounds have been developed for $(R_{0},R_{1})$ configurations (of Definition 2.4), enabling our analysis to cover a much broader range of configurations. To relate these to estimates for (25), we require the following continuity bound for the single layer potential.

Lemma 3.4**.**

For a domain $D$ with bounded Lipschitz boundary $\partial D$, given $k_{0}>0$ the following bound on the single layer potential (5) holds

[TABLE]

Proof 3.5**.**

It is straightforward to show (see, e.g., Hewett et al. (2013, Lemma 4.1)) that

[TABLE]

We shall exploit the Lipschitz property of $\partial D$, by defining a finite set of Lipschitz graphs which describe its geometry, and bounding the right-hand side of (27) in terms of the coordinates describing these graphs. Let $\{W_{j}\}$, $j=1,\ldots,N$, be a finite open cover of $\partial D$ as in the definition of a Lipschitz domain (see, e.g., McLean (2000, 3.28)). Assume without loss of generality that each $W_{j}\cap\partial D$ is connected. Each $W_{j}\cap\partial D$ is part of the graph of a Lipschitz real function $\ell_{j}$ in rotated Cartesian coordinates, which we denote $(x_{j},y_{j})$. The boundary $\partial D$ can thus be decomposed into $N_{D}$ arcs $\alpha_{j}$ (with disjoint relative interiors) that are the graph of $\ell_{j}:[a_{j},b_{j}]\to\mathbb{R}$, i.e. $\alpha_{j}=\{(x_{j},y_{j})\in\mathbb{R}^{2}:\;a_{j}\leq x_{j}\leq b_{j},\;y_{j}=\ell_{j}(x_{j})\}\subset W_{j}\cap\partial D$, and $\partial D=\bigcup_{j=1}^{N_{D}}\alpha_{j}$. Denote by $C_{\ell}$ a constant which bounds above the Lipschitz constant of every Lipschitz graph function $\ell_{j}$. Fix any $\mathbf{p}\in\mathbb{R}^{2}$. For each $j=1,\ldots,N_{D}$ denote by $(p_{x,j},p_{y,j})$ the coordinates of $\mathbf{p}$ in the $(x_{j},y_{j})$ coordinate system. We have $\max\{|a_{j}-p_{x,j}|,|b_{j}-p_{x,j}|\}\leq\max_{\mathbf{q}\in\partial D}|\mathbf{p}-\mathbf{q}|$.

Now we have established the necessary notation, we decompose the integral in the $L^{2}(\partial D)$ norm on the right-hand side of (27) into the regions contained within the open sets $W_{j}$, each with its own Lipschitz graph $\alpha_{j}$:

[TABLE]

Now we may appeal to the monotonicity of $|H^{(1)}_{0}|$, and bound the variation of the mapping to the Lipschitz graph $\ell^{\prime}$ by the constant $C_{\ell}$ to obtain

[TABLE]

where we have changed integration variables to simplify the integrand in the second step. Since $b_{j}-a_{j}\leq R_{D}:=\mathop{\rm diam}(\partial D)$, we can bound further

[TABLE]

We have split the integrals in order to bound the Hankel function, using $|H^{(1)}_{0}(z)|\leq\hat{c}(1+|\log|z||)$ if $0<|z|\leq 1$ with (28), and $|H^{(1)}_{0}(z)|\leq\hat{c}|z|^{-1/2}$ if $|z|>1$, by e.g. Hewett et al. (2013, p. 638) (with value $\hat{c}\approx 2.09$) with (29). The integral (28) is therefore bounded above by

[TABLE]

where we have used $\int^{t}(1-\log s)^{2}\operatorname{d}\!s=t(\log^{2}t-4\log t+5)+$constant in the final step. The second integral (29) is maximised either when (i) $k(a_{j}-p_{x,j})=1$ or when (ii) $k(a_{j}-p_{x,j})=-kR_{D}/2$. In case (i), the integral is bounded above by

[TABLE]

and in case (ii) it is bounded above by

[TABLE]

so in either case, (29) is bounded above by $2\hat{c}^{2}\log(1+kR_{D})$. Combining this with (30) yields

[TABLE]

This gives the explicit form of the simplified estimate in our claim, proving the assertion.

Using this result, we can say more about the $k$-dependence of $u_{\max}(k)$, for a large class of multiple scattering configurations of interest.

Theorem 3.6**.**

Suppose that $u$ satisfies the BVP (2)–(4), with plane wave incidence (1). Then given $k_{0}{>0}$ independent of $k$, the following bounds hold:

1. (i)

If \Upsilon={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\cup{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} is a non-trapping polygon (in the sense of Definition 2.3),

[TABLE] 2. (ii)

Otherwise, if \Upsilon={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\cup{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} is an $(R_{0},R_{1})$ domain (in the sense of Definition 2.4),

[TABLE]

where $u_{\max}(k)$ is as in (25).

Proof 3.7**.**

We write the BVP (2)–(4) for the scattered field $u^{s}$, with Dirichlet data $u^{s}=-u^{i}$ on the boundary $\partial D$, in terms of the Dirichlet-to-Neumann (DtN) map $P_{\text{DtN}}$ (see, e.g., Chandler-Wilde et al. (2012b, §2.7)) as ${\partial u^{s}}/{\partial\mathbf{n}}=-P_{\text{DtN}}\tau_{+}u^{i}$, where $\tau_{+}$ denotes the exterior Dirichlet trace. The representation (6) gives

[TABLE]

This, together with $|{\partial u^{i}}/{\partial\mathbf{n}}|\leq k|u^{i}|$ (which follows immediately from (1)), enables us to bound $u^{s}$ as

[TABLE]

where $\|\cdot\|_{H^{1}_{k}(\partial D)}$ denotes the the $k$-weighted norm of the Sobolev space $H^{1}(\partial D)$

[TABLE]

and $\nabla_{S}$ denotes the surface gradient operator on $\partial D$ (defined in (71)). By the triangle inequality we have $u_{\max}(k)\leq\|u^{i}\|_{L^{\infty}(D)}+\|u^{s}\|_{L^{\infty}(D)}$, and from Lemma 3.4 we can bound $\|{S}_{k}\|_{L^{2}(\partial D)\shortrightarrow L^{\infty}(D)}$. Hence we may write, for $k\geq k_{0}$,

[TABLE]

For the DtN maps, we may use Baskin et al. (2016, Theorem 1.4) for the non-trapping polygon case (i) $\|P_{\text{DtN}}\|_{H^{1}_{k}(\partial D)\shortrightarrow L^{2}(\partial D)}\lesssim 1$, whilst the $(R_{0},R_{1})$ obstacle case (ii) $\|P_{\text{DtN}}\|_{H^{1}_{k}(\partial D)\shortrightarrow L^{2}(\partial D)}\lesssim k^{2}$ follows by Chandler-Wilde et al. (2020, Theorem 1.8). It remains to bound the incident field $u^{i}$ at the boundary and in the domain. For plane wave incidence, it follows by the definitions (32) and (1) that $\|u^{i}\|_{H^{1}_{k}(\partial D)}\leq 2\sqrt{|\Gamma|}k$ and $\|u^{i}\|_{L^{\infty}(D)}=1$. The result follows by combining these bounds on $u^{i}$ with the components of (33).

Theorem 3.6 is a generalisation of Hewett et al. (2013, Theorem 4.3), which bounds $u_{\max}(k)$ for star-shaped polygons. Although more general, Theorem 3.6 differs from Hewett et al. (2013, Theorem 4.3) in that it is not fully explicit in terms of the geometric parameters of $\Upsilon$. We do not expect such a bound to hold for the most general configurations and incident fields, since it was shown in Betcke et al. (2011, Theorem 2.8) that there exist multiple obstacle configurations for which $\|{\mathcal{A}}_{k,\eta}^{-1}\|_{L^{2}(\partial D)\shortrightarrow L^{2}(\partial D)}$ is bounded below by a term which grows exponentially with $k$, in which case $u_{\max}(k)$ would grow similarly. In particular though, Theorem 3.6(i) is immediately applicable to the case of polygons which are non-convex, non-star-shaped and non-trapping, considered in Chandler-Wilde et al. (2015) (see Definition 3.1 therein), for which the stronger result $u_{\max}(k)=\mathcal{O}(1)$ for $k\to\infty$ was conjectured, in the (then) absence of any available algebraic bounds. The bound of Theorem 3.6 is sufficient to guarantee algebraic growth of $u_{\max}(k)$ in $k$, and therefore exponential convergence of HNA-BEM for such polygons.

The following assumption generalises Theorem 3.6 to all configurations of interest.

Assumption 3.1**.**

For the solution $u$ of the BVP (2)–(4), we assume that there exist $\beta\geq 0$, $k_{0}>0$ and ${C_{u}}>0$, independent of $k$, such that

[TABLE]

that is $u_{\max}(k)$ of (25) has at most algebraic dependence on the wavenumber $k$.

Clearly Assumption 3.1 holds for configurations satisfying the conditions of Theorem 3.6 (see Remark 5.5 for more details).

3.3 Analyticity and bounds for the envelope functions $v^{\pm}_{j}$

Additional notation is required for the estimates that follow. Denote by $\Omega_{j}$ the exterior angle at the corner $\mathbf{P}_{j}$ of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} (see figure 3 for an illustrative example). Since {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} is a convex polygon, $\Omega_{j}\in(\pi,2\pi)$ for all $j=1,\ldots,{\mathcal{N}_{\Gamma}}$. Let $c_{*}>0$ be a constant such that $kL_{j}\geq c_{*}$ for all $j=1,\ldots,{\mathcal{N}_{\Gamma}}$ (e.g. $c_{*}=\min_{j=1,\ldots,{\mathcal{N}_{\Gamma}}}\{kL_{j}\}$).

We now aim to show, as in Hewett et al. (2013) where only one (convex polygonal) scatterer {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} is present, that the functions $v_{j}^{\pm}$ are complex-analytic, and moreover that they can be approximated much more efficiently than $\partial u/\partial\mathbf{n}|_{\Gamma}$. We update this to the multiple scattering configuration by adapting the intermediate results of Hewett et al. (2013, §3). We first consider the solution behaviour near the corners.

Lemma 3.8** (Solution behaviour near the corners).**

Suppose that $u$ satisfies the BVP (2)–(4) and $\mathbf{x}\in D$ satisfies $r:=|\mathbf{x}-\mathbf{P}_{j}|\in(0,1/k]$, and $r<\mathop{\rm dist}(\mathbf{P}_{j},\gamma)$. Then there exists a constant $C>0$, depending only on $\partial D$ and $c_{*}$, such that (with $u_{\max}(k)$ as in (25)),

[TABLE]

Proof 3.9**.**

Follows identical arguments to Hewett et al. (2013, Lemma 3.5), with the slight modification to the definition $R_{j}:=\min\{L_{j-1},L_{j},\pi/(2k),\mathop{\rm dist}(\mathbf{P}_{j},\gamma)\}$, which ensures only areas close to the corner $\mathbf{P}_{j}$ inside $D$ are considered.

Now we may bound the singular behaviour of the diffracted envelopes $v_{j}^{\pm}$, which will enable us to choose a suitable approximation space for the numerical method.

Theorem 3.10**.**

Suppose that $u$ is a solution of the BVP (2)–(4), and that $c_{r}\in(0,1]$ is chosen such that $\mathop{\rm dist}(\{\mathbf{P}_{j}:{j=1,\ldots,{\mathcal{N}_{\Gamma}}}\},\gamma)>c_{r}/k$. Then the diffracted wave envelope components $v_{j}^{\pm}$ for $j=1,\ldots,{\mathcal{N}_{\Gamma}}$, of the boundary representation (3.1), are analytic in the right complex half-plane $\mathop{\rm Re}[s]>0$, where they satisfy the bounds

[TABLE]

where $\delta^{+}_{j},\delta^{-}_{j}\in(0,1/2)$ are given by $\delta^{+}_{j}:=1-\pi/\Omega_{j}$ and $\delta^{-}_{j}:=1-\pi/\Omega_{j+1}$. The constant $C^{+}_{j}$ depends only on $c_{*}$, $c_{r}$ and $\Omega_{j}$, whilst the constant $C^{-}_{j}$ depends only on $c_{*}$, $c_{r}$ and $\Omega_{j+1}$.

Proof 3.11**.**

The analyticity of the functions $v^{\pm}_{j}(s)$ in $\mathop{\rm Re}[s]>0$ follows from their definition (20)–(21) and the analyticity of $\mu(s)$ in the same set, which is shown in Hewett et al. (2013, Lemma 3.4). The estimate of $|v_{j}^{\pm}(s)|$ for $|s|>1/k$ follows as in the proof of Hewett et al. (2013, Theorem 3.2). Here we show for $v^{+}_{j}$, the proof for $v^{-}_{j}$ follows similar arguments. For $|s|\leq 1/k$, the definition (20)–(21) of $v_{j}^{+}$ gives

[TABLE]

Since $c_{r}\leq 1$ and thanks to Lemma 3.8, the first integral is bounded as in the proof of Hewett et al. (2013, Theorem 3.2), leading to the term $u_{\max}(k)k|ks|^{-\delta^{\pm}_{j}}$ in the assertion. Using the bound on $\mu$ from Hewett et al. (2013, Lemma 3.4), we control the second integral as

[TABLE]

The bound in the assertion follows by noting that $k|s|+c_{r}<2$.

The constant $c_{r}$ is small when the scatterers are close together, relative to the wavelength of the problem. Thus the terms containing $c_{r}$ in the bound of Theorem 3.10 control the effect of the separation between {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} on the singular behaviour of $v^{\pm}_{j}$. However, the method we present is designed for high-frequency problems, and to maintain $c_{r}=O(1)$ as $k$ increases, the separation of the scatterers is allowed to decrease inversely proportional to $k$. Hence, for the configurations that we consider of practical interest in the high-frequency regime, the condition (3.12) in the following corollary will hold.

Corollary 3.12**.**

Suppose that the conditions of Lemma 3.8 hold, with the additional constraint that the separation condition

[TABLE]

is satisfied. It then follows that the first bound of Theorem 3.10 can be simplified to

[TABLE]

Proof 3.13**.**

If the separation condition (34) holds, we can choose $c_{r}=1$ in Theorem 3.10, from which $(k|s|+c_{r})^{-1}\leq 1$. The term $k(k|s|+c_{r})^{-1}$ is therefore dominated by the term $k|ks|^{\delta^{\pm}_{j}}$ for $0<|s|\leq 1/k$.

The separation condition (34) aligns the bounds of Theorem 3.10 with the well-studied single scattering HNA configurations of Hewett et al. (2013, Theorem 5.2). Hence, all best approximation results for the single scattering case may be applied to the approximation on $\Gamma$ in the multiple scattering problems we consider here.

Remark 3.14**.**

The result of Theorem 3.10 may be extended to an incident wave of source-type, for example the point source emanating from $\mathbf{s}\in D$, $u^{i}(\mathbf{x})=H^{(1)}_{0}(k|\mathbf{x}-\mathbf{s}|)$. This requires that the position of the source point $\mathbf{s}$ is separated by a distance of at least $1/k$ from {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} (similar to the separation condition (34)), see Gibbs (2017, §3.2) for details.

4 $hp$ approximation space

We will combine two approximation spaces: the HNA-BEM space on $\Gamma$ and a standard $hp$-BEM space on $\gamma$. Hereafter, using the parametrisation of the boundaries $\Gamma$ and $\gamma$, we identify $L^{2}(\Gamma_{j})$ with $L^{2}(0,L_{j})$, and $L^{2}(\gamma)$ with $L^{2}(0,L_{\gamma})$.

4.1 HNA-BEM approximation on $\Gamma$

As in previous HNA methods, on $\Gamma$ we approximate only the diffracted waves

[TABLE]

where $v_{j}^{\pm}$ are as in (20)–(21), and broadly speaking this is done using basis elements of the form

[TABLE]

where $P_{j}^{\pm}$ are piecewise polynomials on a graded mesh. There are two well-studied classes of $hp$ approximation space we may use to do this. Both spaces consist of piecewise polynomials multiplied by oscillatory functions oscillating in both directions along the surface of $\Gamma$, and both spaces are constructed on meshes graded towards the singularities at the corners of $\Gamma$. We briefly describe these approximation spaces here:

1. (i)

The overlapping-mesh space, used in original HNA methods for single scatterers, this discrete space is the sum of two subspaces, each constructed on a separate mesh graded in opposite directions. On $\Gamma_{j}$, the subspace on the mesh graded towards ${\widetilde{L}}_{j-1}$ is used to approximate $v_{j}^{+}(s-{\widetilde{L}}_{j-1})\mathrm{e}^{\mathrm{i}ks}$ and the subspace on the mesh graded towards ${\widetilde{L}}_{j}$ is used to approximate $v_{j}^{-}({\widetilde{L}}_{j}-s)\mathrm{e}^{-\mathrm{i}ks}$. Details can be found in Hewett et al. (2013, §5). 2. (ii)

The single-mesh space, constructed on a single mesh graded towards both edges. This space can easily be implemented by adapting a standard BEM code, as the mesh is of a more standard type. However, care must be taken close to the corners of $\Gamma$: certain elements must be removed from the approximation space to ensure the discrete system does not become too ill-conditioned. We will define this space shortly.

A range of numerical experiments comparing both approximation spaces for collocation HNA-BEM can be found in Parolin (2015). For either choice of mesh, we denote by ${n_{j}}$ the number of grading layers and by $p_{j}$ the maximum polynomial degree on the $j$th side (in terms of the notation of Chandler-Wilde & Langdon (2007) and Hewett et al. (2013), we choose $p_{j}=p^{+}_{j}=p^{-}_{j}$, ${n_{j}}={n_{j}}^{+}={n_{j}}_{-}$ for simplicity). We denote by $\sigma>0$ the grading parameter, so that the smallest mesh element of $\Gamma_{j}$ (touching the corners of $\Gamma_{j}$) has length $L_{j}\sigma^{n_{j}}$.

The single-mesh space has been described in the theses Gibbs (2017); Parolin (2015) and is used for the numerical experiments in §6; we define it here for convenience.

Definition 4.1**.**

Given ${L}>0$, $n\in\mathbb{N}$ and a grading parameter $\sigma\in(0,1/2)$, we denote by ${\mathcal{M}}_{n}(0,{L})=\{x_{0},\ldots,x_{2n+1}\}$ the symmetric geometrically graded mesh on $[0,{L}]$ with $n$ layers in each direction, whose $2n+2$ meshpoints $x_{i}$ are defined by

[TABLE]

For a vector $\mathbf{p}=(p_{1},\ldots,p_{n+1})\in(\mathbb{N}_{0})^{n+1}$ we denote by ${\mathcal{P}}_{\mathbf{p},n}$ the space of piecewise polynomials on ${\mathcal{M}}_{n}(0,{L})$ with degree vector $\mathbf{p}$, i.e.

[TABLE]

We first define two spaces for each side $\Gamma_{j}$, $j=1,\ldots,{\mathcal{N}_{\Gamma}}$, using ${n_{j}}\in\mathbb{N}$ to determine the degree of mesh grading and the vectors $\mathbf{p}_{j}$ to determine the polynomial degree on each mesh element:

[TABLE]

As is explained in Remark 4.2, to avoid ill-conditioning of the discrete system we must remove certain basis functions supported on the elements within a given distance from the corners:

[TABLE]

where

[TABLE]

and $\alpha_{j}$ is a parameter chosen such that $0<\alpha_{j}<L_{j}k/(4\pi)$, bounded independently of $k$ and $\mathbf{p}_{j}$, used to fine-tune the space. Put simply, there are two basis functions on (large) elements sufficiently far from the corners, and one basis element on (small) elements close to the corners. The parameter $\alpha_{j}$ determines what is meant by sufficiently far. Hence the single-mesh approximation space with dimension ${N_{\Gamma}}$ is defined as

[TABLE]

Remark 4.2** (Why basis elements of the single-mesh space are removed).**

Since the mesh is strongly graded to approximate the singularities of $v^{\pm}_{j}$, some of its elements are much smaller than the wavelength of the problem, thus on these elements $\mathrm{e}^{\pm\mathrm{i}ks}$ are roughly constant and the functions of ${V}^{+}_{j}$ supported on these elements are numerically indistinguishable from those on ${V}^{-}_{j}$, leading to an ill-conditioned discrete system of Galerkin equations set in ${V}_{j}^{+}\cup{V}_{j}^{-}$. To avoid this, in these elements we maintain only one of these two contributions. Intuitively, $\alpha_{j}$ can be thought of as the value such that in all mesh elements with distance from one of the segment endpoints smaller than $\alpha_{j}$, the space ${\widetilde{V}}_{j}$ supports polynomials multiplied with only one of the waves $\mathrm{e}^{\pm\mathrm{i}ks}$. As the parameter $\alpha_{j}$ increases, fewer degrees of freedom are used and the conditioning of the discrete system is improved, but the accuracy of the method is reduced, hence care must be taken when selecting $\alpha_{j}$.

In much of what follows, the choice of single- or overlapping-mesh HNA space is irrelevant, hence we shall use $V_{{N_{\Gamma}}}^{\mathop{\rm HNA}}(\Gamma)$ to denote either, but will make clear the cases for which the choice is significant. For the overlapping-mesh space, best approximation estimates were derived in Hewett et al. (2013, Theorem 5.4). The following result from Gibbs (2017, Corollary 2.11) compares the best approximation of the single-mesh and overlapping-mesh spaces, on $\Gamma$.

Theorem 4.3**.**

Suppose that the obstacles {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} are sufficiently far apart so that the separation condition (34) holds. Let $V_{{N_{\Gamma}}}^{\mathop{\rm HNA}}(\Gamma)$ be an HNA space as above, $c_{j}>0$ be such that the polynomial degrees $p_{j}$ and the numbers of layers $n_{j}$ satisfy

[TABLE]

and denote $p_{\Gamma}:=\min_{j}\{p_{j}\}$. Then we have the following best approximation estimate for the diffracted wave $v_{\Gamma}$ (of (35)):

[TABLE]

where $C_{\Gamma}$ is a constant independent of $k$ and

[TABLE]

with $I_{*}$ and $\tau_{\Gamma}$ independent of $n_{j},p_{j},k$ (both are defined precisely in Gibbs (2017, Corollary 2.11)), $\delta_{*}:=\min_{j,\pm}\{\delta^{\pm}_{j}\}$ (with $\delta_{j}^{\pm}$ as in Theorem 3.10), whilst $L_{*}:=\max_{j}L_{j}$ the length of the longest side of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}. For the single-mesh space, it follows that $C_{\Gamma}=\max_{j}\{C_{j}\}$ for $C_{j}$ of Gibbs (2017, Theorem 2.9). For the overlapping-mesh space, $C_{\Gamma}$ is equal to the constant $C_{4}$ of Hewett et al. (2013, Theorem 5.5).

Theorem 4.3 shows that we obtain exponential convergence of the best approximation to $v_{\pm}$ with respect to $p_{\Gamma}$, which controls both polynomial degree and mesh grading (via (36)), across all wavenumbers $k$. To maintain accuracy as $k$ increases one needs to increase $p_{\Gamma}$ in proportion to $\log k$, and hence the total number of degrees of freedom (which is proportional to $p_{\Gamma}^{2}$) in proportion to $\log^{2}k$.

Remark 4.4**.**

It is shown in Hewett et al. (2011, Theorem A.3) for the overlapping-mesh HNA space that it is possible to reduce the number of degrees of freedom on $\Gamma$, whilst maintaining exponential convergence, by reducing the polynomial degree in the smaller mesh elements, as is standard in $hp$ schemes. For example, given a polynomial degree $p_{j}>1$, we can define for each side $\Gamma_{j}$, $j=1,\ldots,{\mathcal{N}_{\Gamma}}$, a degree vector $\mathbf{p}_{j}$ by

[TABLE]

where $n_{j}$ is as in Definition 4.1 of the single-mesh space. This may be applied to either the single or overlapping mesh, and results in a linear reduction of polynomial degree on mesh elements closer to the corners of $\Gamma_{j}$. Numerical experiments in §6 suggest that exponential convergence is maintained for the single-mesh HNA space if the degrees of freedom are reduced in this way, although we do not prove this here.

4.2 Standard $hp$-BEM approximation on $\gamma$

If Assumption 3.1 holds, as is the case in the configurations of Theorem 3.6, it follows from Theorem 4.3 that it is sufficient for the number of DOFs in $V_{{N_{\Gamma}}}^{\mathop{\rm HNA}}(\Gamma)$ to grow logarithmically with $k$, to accurately approximate $v_{\pm}$. However, this tells us nothing about the DOFs required on $\gamma$. To account for the contribution from $\gamma$, we parametrise $\mathbf{x}_{\gamma}:[0,L_{\gamma}]\rightarrow\gamma$ and construct an appropriate (depending on the geometry of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}) ${N_{\gamma}}$-dimensional approximation space $V_{{N_{\gamma}}}^{{hp}}(\gamma)\subset L^{2}(0,L_{\gamma})$ for

[TABLE]

While a representation analogous to (3.1) holds on $\gamma$ when {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} is a convex polygon, this approach is not suitable for the present multiple scattering approximation. If such a representation were used on multiple polygons, the system to solve would need to be written as a Neumann series and solved iteratively. This alternative approach is outlined briefly in Gibbs (2017, §4.4.1). Instead we approximate the full solution $v_{\gamma}$, rather than any of its individual components as listed in (3.1). An advantage of the approach in this paper is that the only restriction imposed on $\gamma$ is that it must be Lipschitz and piecewise analytic. The disadvantage is that the number of DOFs required to approximate the solution on $\gamma$ has to increase with frequency to maintain accuracy, as is typical of standard $hp$-schemes. For all $k$ such that $|\gamma|$ is small compared with the wavelength $2\pi/k$, one would not expect this increase in DOFs to be significant. Here, we take $V_{{N_{\gamma}}}^{{hp}}(\gamma)$ to be a standard $hp$-BEM approximation space consisting of piecewise polynomials to approximate $v_{\gamma}$, with mesh and degree vector dependent on the geometry of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}.

We now aim to bound the approximation of the solution on $\gamma$, in terms of key parameters, for the case where $\gamma$ is analytic. This will enable us to quantify the $k$-dependence of our method, which we expect to be mild when $|\gamma|$ is small compared with the wavelength. A range of tools were developed in Löhndorf & Melenk (2011) for $hp$-BEM approximations for problems of scattering by analytic surfaces, provided bounds on ${\mathcal{A}}_{k,\eta}^{-1}$ are available. For this, we are able to use recently developed theory of $(R_{0},R_{1})$ configurations (of Definition 2.4) for which we have from Chandler-Wilde et al. (2020, (1.28)): if $\eta=O(k)$, then given $k_{0}>0$,

[TABLE]

In the class of problems we consider, the total boundary $\Gamma\cup\gamma$ is not analytic, because $\Gamma$ is the boundary of a polygon. Therefore we could not apply the theory of Löhndorf & Melenk (2011) to a standard $hp$ approximation on $\Gamma\cup\gamma$. However, in our method the standard $hp$ approximation is only on $\gamma$, which in this section we will restrict to be analytic; Theorem 4.3 provides a best approximation estimate for the HNA space on the polygon $\Gamma$. As we shall see, this is sufficient to get a best approximation estimate for $v_{\gamma}$ in the standard $hp$ space $V_{{N_{\gamma}}}^{{hp}}(\gamma)$. The main idea is to consider an equivalent problem of scattering by (only) the obstacle {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}, with the contribution from {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} absorbed into the incident field. We can rewrite the representation (6)

[TABLE]

separating the contribution from the convex polygon $\Gamma$. To construct an equivalent problem, we consider the additional component of the incident field to be the contribution from $\Gamma$:

[TABLE]

where $T_{\gamma}$ is a tubular neighbourhood of $\gamma$, i.e. for some $\epsilon>0$ we have

[TABLE]

with $\epsilon$ chosen such that $\mathop{\rm dist}(T_{\gamma},\Gamma)>0$. Our equivalent problem is therefore scattering of $u^{i}+u^{i}_{\Gamma}$ by {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}, in $T_{\gamma}$. It is straightforward to see that the solution to this equivalent problem is the same as the solution to the BVP (2)–(4) (restricted to $T_{\gamma}$). To use the $hp$ theory developed in Löhndorf & Melenk (2011), we must show that the solution to our scattering problem is in the space of Löhndorf & Melenk (2011, Definition 1.1):

[TABLE]

for some $\xi$ independent of $k$, $h$, $p$ and

[TABLE]

A prerequisite for $u\in\mathcal{U}(\psi,\xi,T_{\gamma}\setminus\gamma)$ is that the incident field to our equivalent problem $u^{i}+u^{i}_{\Gamma}$ is also in $\mathcal{U}(\psi,\xi,T_{\gamma}\setminus\gamma)$, possibly for different parameters $\psi$ and $\xi$.

Lemma 4.5**.**

If {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\cup{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} is an $(R_{0},R_{1})$ configuration, then

[TABLE]

where $\psi(k):=Ck^{7/2}\log^{1/2}(k\mathop{\rm diam}(\Gamma)+1)$ with $C>0$ a constant independent of $k$.

Proof 4.6**.**

Throughout the proof we let $C$ denote an arbitrary constant independent of $k$ and $n$. It follows from standard mapping properties of the single-layer operator (e.g., Chandler-Wilde et al. (2012b, Theorem 2.15(i))) that $u^{i}_{\Gamma}\in H^{1}(\omega)$, where $\omega$ is a bounded open subset of $\mathbb{R}^{2}$ containing {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\cup{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}. We may therefore bound using Melenk (2012, Theorem B.6), choosing zero forcing term to obtain

[TABLE]

given $k_{0}>0$, where $\omega$ is a bounded open set compactly containing $T_{\gamma}$ and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}. From (41), we see that the norm is the sum of $n+1$ terms, hence

[TABLE]

given $k_{0}>0$, which follows by combining with (42) and $(n+1)!\leq(n+1)^{n}$. We now bound $u^{i}_{\Gamma}$ in terms of known quantities,

[TABLE]

We may bound these norms using Lemma 26, (38) and (8) (choosing $\eta=O(k)$) to obtain

[TABLE]

Finally, we can combine the bound (45) with (44) to obtain

[TABLE]

proving the assertion.

Now we have shown sufficient conditions on the growth of the derivatives of $u^{i}_{\Gamma}$, we are ready to obtain best approximation estimates on $\gamma$.

Proposition 4.7**.**

Suppose \Upsilon={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\cup{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} is an $(R_{0},R_{1})$ configuration (in the sense of Definition 2.4) and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} has an analytic boundary $\gamma$. If $V_{{N_{\gamma}}}^{{hp}}(\gamma)$ is constructed on a quasi-uniform mesh (in the sense of Löhndorf & Melenk (2011, §1)) with $kh/p_{\gamma}\lesssim 1$, where $h$ and $p_{\gamma}$ denote maximum mesh width and polynomial degree respectively, then given positive constants $k_{0}$, $\zeta$ independent of $k,p_{\gamma}$ and $h$ we have the following best approximation estimate:

[TABLE]

where

[TABLE]

with $C>0$ a constant independent of $k,p_{\gamma}$ and $h$.

Proof 4.8**.**

By Lemma 4.5 we have that $u^{i}_{\Gamma}\in\mathcal{U}(\psi,1,T_{\gamma}\setminus\gamma)$, and it is straightforward to see that $u^{i}\in\mathcal{U}(1,1,T_{\gamma}\setminus\gamma)$. Choosing

[TABLE]

with $g_{1}=g_{2}=0$ otherwise, we have that $g_{1},g_{2}\in\mathcal{U}(\psi,1,T_{\gamma}\setminus\gamma)$. Noting again (38), we may appeal to Löhndorf & Melenk (2011, Lemma 2.6) to deduce that the solution $u$ of the BVP (2)–(4) (which is the same as the solution to the equivalent problem of scattering by $u^{i}+u^{i}_{\Gamma}$), is in the space $\mathcal{U}(\psi^{*},1,T_{\gamma}\setminus\gamma)$, for all $k\geq k_{0}$ given $k_{0}>0$, where

[TABLE]

Hence, the best approximation estimate of Löhndorf & Melenk (2011, Lemma 3.16) may be applied to

[TABLE]

(noting Löhndorf & Melenk (2011, Definition 3.3)), yielding the best approximation result after rescaling by $C_{\gamma}(k)$.

We do not expect the above result to be sharp, however to the best knowledge of the authors, it is the only $hp$-BEM estimate currently available for such a configuration. We now generalise Proposition 4.7 in the form of an assumption, which states that we observe exponential convergence to the solution $\partial u/\partial\mathbf{n}$ in $V_{{N_{\gamma}}}^{{hp}}(\gamma)$. It follows immediately from Proposition 4.7 that this assumption holds for analytic $\gamma$, under appropriate conditions. For the case of polygonal $\gamma$, the numerical experiments of §6 suggest the assumption also holds, provided that we fix ${N_{\gamma}}=O(k)$.

Assumption 4.1**.**

Denoting by $v_{\gamma}$ the restriction to $\gamma$ of the solution of the BIE (7), we assume that the sequence of approximation spaces

[TABLE]

is such that

[TABLE]

where the positive constants $C_{\gamma}(k),\tau_{\gamma}(k)$ may depend on $k$, and $p_{\gamma}$ is the polynomial degree of the space $V_{{N_{\gamma}}}^{{hp}}(\gamma)$.

4.3 Combined approximation space on $\Gamma\cup\gamma$

The approximation space is based on the representation of the Neumann trace

[TABLE]

where $v_{\Gamma}$ and $v_{\gamma}$ are the unknowns that we solve for using the approximation spaces of §4.1 and §4.2, whilst $\Psi$ denotes the Physical Optics Approximation (17) and $\mathcal{G}_{\gamma\shortrightarrow\Gamma}$ denotes the Interaction Operator (22). Hence the approximation lies in the space

[TABLE]

where the total number of degrees of freedom is $N={N_{\Gamma}}+{N_{\gamma}}$. For problems of one large polygon and one (or many) small polygon(s), the single-mesh HNA space $V_{{N_{\Gamma}}}^{\mathop{\rm HNA}}(\Gamma)$ is particularly practical, as only a small modification is required to implement both this and a standard $hp$-BEM space on $V_{{N_{\gamma}}}^{{hp}}(\gamma)$.

The following notation will be used to describe the problem in block operator form.

Definition 4.9** (Operator restriction).**

For the operator ${\mathcal{A}}_{k,\eta}:L^{2}(\partial D)\to L^{2}(\partial D)$ (of Definition 2.1) and relatively open $X,Y\subset\partial D$, we define the operator $\mathcal{A}_{Y\shortrightarrow X}:L^{2}(Y)\rightarrow L^{2}(X)$ by

[TABLE]

where $\mathcal{Q}_{Y}:L^{2}(Y)\to L^{2}(\partial D)$ is the zero-extension operator, such that $(\mathcal{Q}_{Y}\varphi)|_{Y}=\varphi$ and $(\mathcal{Q}_{Y}\varphi)|_{\partial D\setminus Y}=0$. For the case of the identity operator $I_{X\shortrightarrow X}:L^{2}(X)\rightarrow L^{2}(X)$ we simplify the notation by writing $I_{X}$.

Inserting (47) into the BIE (7), we can write the problem to solve in block form: Find $v\in L^{2}(\Gamma)\times L^{2}(\gamma)$ such that

[TABLE]

Stated in a variational form equivalent to (9), our problem is as follows: Find $v\in L^{2}(\Gamma\cup\gamma)$ such that

[TABLE]

for all $w\in L^{2}(\Gamma\cup\gamma)$. This problem forms the basis of our Galerkin method.

5 Galerkin method

In this section, we derive error bounds for the approximation of equation (51) by the Galerkin method on the discrete space $V_{N}^{\mathop{\rm HNA}^{*}}(\Gamma,\gamma)$ (defined in (48)). Under certain assumptions, we will show that exponential convergence is achieved. We intend to approximate the unknown components of the solution on $\Gamma$ and $\gamma$, that is

[TABLE]

where $v$ is the solution to (51). Recall (from §4.3) that we use an HNA approximation space $V_{{N_{\Gamma}}}^{\mathop{\rm HNA}}(\Gamma)$ (single- or overlapping-mesh) on $\Gamma$, with a standard $hp$-approximation space $V_{{N_{\gamma}}}^{{hp}}(\gamma)$ on $\gamma$. The discrete problem to solve is: find $v_{N}\in V_{N}^{\mathop{\rm HNA}^{*}}(\Gamma,\gamma)$ such that

[TABLE]

for all {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}({w^{N}_{\Gamma}},{w^{N}_{\gamma}})}\in V_{N}^{\mathop{\rm HNA}^{*}}(\Gamma,\gamma). To implement the Galerkin method, we choose suitable bases $\Lambda_{\Gamma}$ and $\Lambda_{\gamma}$, with

[TABLE]

To determine $v_{N}$ we seek ${\mathbf{a}}\in\mathbb{C}^{N}$ which solves the block matrix system $B\mathbf{a}=\mathbf{b}$, where

[TABLE]

and

[TABLE]

For further details on implementation, see Remark 6.1.

For the remainder of the section, we present approximation estimates of quantities of practical interest. We assume that as $N$ increases, so do ${N_{\Gamma}}$ and ${N_{\gamma}}$, such that the following convergence conditions hold:

[TABLE]

for $i=1,\ldots{\mathcal{N}_{\gamma}}$, where $w_{N}:=(w_{\Gamma}^{N},w_{\gamma_{1}}^{N},\ldots,w_{\gamma_{\mathcal{N}_{\gamma}}}^{N})$. If $V_{{N_{\gamma}}}^{{hp}}(\gamma)$ is a standard $hp$-BEM space then (64) holds. It follows by identical arguments to Chandler-Wilde & Langdon (2007, Theorem 5.1) and the definition of the single- and overlapping-mesh spaces (see §4) that Condition (63) holds. We present a lemma concerning the stability of the system (59)–(60).

Lemma 5.1** (Stability of discrete system).**

Suppose the convergence conditions (63)–(64) hold. Then there exist positive constants $C_{q}(k)$ and $N_{0}(k)$ such that for $N\geq N_{0}$ the solution $v_{N}$ of (59)–(60) exists. Moreover

[TABLE]

Proof 5.2**.**

First we show that $\mathcal{A}_{\square}$ is a compact perturbation of an operator which is Fredholm of zero index. We have from Chandler-Wilde & Langdon (2007, p. 620) that $\mathcal{A}_{\Gamma\shortrightarrow\Gamma}$ is a compact perturbation of a Fredholm operator (of index zero), and the same arguments can be applied to each $\mathcal{A}_{\gamma_{i}\shortrightarrow\gamma_{i}}$ for $i=1,\ldots{\mathcal{N}_{\gamma}}$. As the kernels of $\mathcal{A}_{{\Gamma\shortrightarrow\gamma}_{i}}$, $\mathcal{A}_{\gamma_{i}\shortrightarrow\Gamma}$ and ${\mathcal{A}}_{\gamma_{i}\shortrightarrow\gamma_{\ell}}$ for $i\neq\ell$ are continuous for $i=1,\ldots,{\mathcal{N}_{\gamma}}$, these operators are also compact, hence $\mathcal{A}_{k,\eta}$ is a compact perturbation of a coercive Fredholm of zero index operator.

Let $\mathcal{P}_{N}$ be the orthogonal projection operator from $L^{2}(\Gamma)\times L^{2}(\gamma)$ onto $V_{N}^{\mathop{\rm HNA}^{*}}(\Gamma,\gamma)$. Given the convergence condition (63), it follows by the density of $C^{\infty}(0,L_{\Gamma})$ in $L^{2}(0,L_{\Gamma})$ for $j=1,\ldots,{\mathcal{N}_{\Gamma}}$ that we have convergence of the best approximation to any $L^{2}(0,L_{\Gamma})$ function in $V_{{N_{\Gamma}}}^{\mathop{\rm HNA}}(\Gamma)$. Similar arguments follow for convergence on $\gamma$, by the convergence condition (64). Then Chandler-Wilde & Langdon (2007, Theorem 5.2) shows the existence of a solution to the discrete problem (59)–(60), for $N$ sufficiently large, via a bound on

[TABLE]

where

[TABLE]

To show that our method converges to the true solution, we proceed as in Chandler-Wilde & Langdon (2007, Theorem 5.3), noting that

[TABLE]

which we combine with (51) to obtain

[TABLE]

Rearranging and adding $v$ to both sides yields

[TABLE]

hence we can bound

[TABLE]

and the bound follows from the definition of $\mathcal{P}_{N}$ and (65).

For our operator ${\mathcal{A}}_{k,\eta}$, there is little that can be said about the constants $C_{q}(k)$ and $N_{0}(k)$ for the scattering configurations considered in this paper. In the appendix we introduce an alternative BIE formulation which is coercive provided that $|\gamma|$ is of the order of one wavelength. For this coercive formulation, $N_{0}(k)=1$ and $C_{q}(k)$ can be made explicit.

Recalling that we are actually approximating the (dimensionless) diffracted waves on $\Gamma$ and the (dimensionless) Neumann trace of the solution on $\gamma$, the full approximation to the Neumann trace follows by inserting $v_{N}$ into (47) and is denoted

[TABLE]

The following theorem can be used to determine the error of the full approximation.

Theorem 5.3**.**

Suppose that

1. (i)

the separation condition (34) holds, 2. (ii)

the convergence conditions (63)–(64) hold, 3. (iii)

Assumption 4.1 (exponential convergence of $V_{{N_{\gamma}}}^{{hp}}(\gamma)$) holds, 4. (iv)

Assumption 3.1 (algebraic growth of the solution of the BVP (2)–(4)) holds.

Then we have the following bound on the error of the approximation (66) to the solution ${\partial u}/{\partial\mathbf{n}}$:

[TABLE]

for $N\geq N_{0}$, where

1. (i)

$N_{0}$* and $C_{q}$ are as in Lemma 5.1,* 2. (ii)

$C_{\mathcal{G}}$* as in Lemma 3.2,* 3. (iii)

$C_{u}$* and $\beta$ are the constants from Assumption 3.1,* 4. (iv)

$p_{\Gamma}$, $J(k)$, $C_{\Gamma}$ and $\tau_{\Gamma}$ are as in Theorem 4.3, 5. (v)

$p_{\gamma}$, $C_{\gamma}$ and $\tau_{\gamma}$ are as in Assumption 4.1.

Proof 5.4**.**

First we focus on the best approximation of ${\partial u}/{\partial\mathbf{n}}$ by an element $w=(w_{\Gamma},w_{\gamma})$ of $V_{N}^{\mathop{\rm HNA}^{*}}(\Gamma,\gamma)$. By the definition (47) we have

[TABLE]

Applying Lemma 5.1 and recalling the definition (48) of $V_{N}^{\mathop{\rm HNA}^{*}}(\Gamma,\gamma)$, we can write

[TABLE]

The assertion follows by combining this inequality with Lemma 3.2, Assumption 3.1, Theorem 4.3 and Assumption 4.1.

For a fixed frequency, Theorem 5.3 suggests that the proposed method is well suited to problems for which {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} is a convex polygon, and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} has a size parameter much smaller than {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}. This is because the number of DOFs required to maintain accuracy in the approximation space on $\Gamma$ grows only logarithmically with $k$. The method will hence be particularly effective if {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} has a size parameter of the order of one wavelength, since in this case the oscillations on $\gamma$ are resolved whilst ${N_{\Gamma}}$ does not need to be large to account for high frequencies due to the (almost) frequency independence of the approximation on $\Gamma$.

Remark 5.5** (Dependencies of parameters of Theorem 5.3).**

In the following situations the bounding constants of Theorem 5.3 can be made either fully explicit, or $k$-explicit.

1. (i)

The terms $C_{\Gamma}$, $\tau_{\Gamma}$ and $J(k)$ are fully explicit given $k$, the geometry of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} and the parameters of $V_{{N_{\Gamma}}}^{\mathop{\rm HNA}}(\Gamma)$. This follows from the separation condition (34). 2. (ii)

In the appendix we present an alternative boundary integral equation which is coercive, under certain geometric restrictions. In such a case $C_{q}(k)$ is known and $N_{0}(k)=1$. 3. (iii)

By Theorem 3.6, if {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\cup{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} is a non-trapping polygon (in the sense of Definition 2.3), then we can choose $\beta=1/2+\epsilon$ for any $\epsilon>0$. 4. (iv)

If $\Upsilon$ is an $(R_{0},R_{1})$ configuration, then by Theorem 3.6 we obtain $\beta=5/2+\epsilon$ for $\epsilon>0$. Furthermore, if $\gamma$ is also analytic and $V_{{N_{\gamma}}}^{{hp}}(\gamma)$ satisfies the conditions of Proposition 4.7 we have $C_{\gamma}(k)=Ck^{8}\log({k\mathop{\rm diam}(\Gamma)+1})$, and $\tau_{\gamma}$ is given by (46).

An approximation $u_{N}$ to the solution $u$ of the BVP (2)–(3) in $D$ is obtained by combining $\nu_{N}$ with the representation formula (6),

[TABLE]

Here the parametrisation $\mathbf{y}_{\Gamma}$ is as in (18) and $\mathbf{y}_{\gamma}$ as in §4.2. Expanding further, we can extend the definition of $\mathcal{G}_{\gamma\shortrightarrow\Gamma}$ to a parametrised form by

[TABLE]

where the indicator function

[TABLE]

is used to ensure the path of integration remains inside the relative upper half-plane $U_{j}$.

Corollary 5.6**.**

Assume conditions (i)–(iv) of Theorem 5.3 hold. Then given $k_{0}~{}>~{}0$, the HNA-BEM approximation to the BVP (2)–(4) satisfies the error bound

[TABLE]

for $N\geq N_{0}$ and $k\geq k_{0}$. The terms in the bound are as in Theorem 5.3.

Proof 5.7**.**

The result follows from the representation (5), the bounds on $\|S_{k}\|_{L^{2}(\partial D)\shortrightarrow L^{\infty}(D)}$ given in Lemma 3.4, Theorem 5.3, and

[TABLE]

A quantity of practical interest is the far-field pattern of the scattered field $u^{s}$, which describes the distribution of energy of the scattered field $u^{s}$ (of a solution to the BVP (2)–(4)) far away from {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\cup{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}. We can represent the asymptotic behaviour of the scattered field (as in Hewett et al. (2013, §6)) by

[TABLE]

where the term $u^{\infty}(\theta)$ denotes the far-field pattern at observation angle $\theta\in[0,2\pi)$, which we can represent via the solution to the BIE (7):

[TABLE]

We may define an approximation $u^{\infty}_{N}$ to the far-field pattern $u^{\infty}$ by inserting $\nu_{N}$ into (68) in place of ${\partial u}/{\partial\mathbf{n}}$.

Corollary 5.8**.**

Under the assumption of Theorem 5.3, the far-field pattern $u^{\infty}_{N}$ computed from the HNA-BEM solution approximates $u^{\infty}$ with the error bound

[TABLE]

The terms in the bound are as in Theorem 5.3.

Proof 5.9**.**

We have

[TABLE]

and the result follows by Theorem 5.3.

6 Numerical results

Here we present numerical results for the solution of the discrete problem (59)–(60). Experiments were run over a range wavenumbers $k\in\{20,40,80,160\}$, incident angles $\mathbf{d}$ and maximal polynomial degrees $p\in\{1,\ldots,8\}$, for three scattering configurations, which we shall refer to as Experiments One, Two and Three. Each configuration consists of an equilateral triangle $\Omega$ with perimeter $L_{\Gamma}=6\pi$, and some small scatterer(s) $\omega$. In Experiment One (§6.1), $\omega$ consists of a single small triangular scatterer with perimeter $L_{\gamma}=3\pi/5$, with the obstacles separated by a fixed distance of $\mathop{\rm dist}(\Gamma,\gamma)=\sqrt{3}\pi/5$, as in Figure 5LABEL:sub@fig:exp1D. In Experiment Two (§6.2) we reduce the distance between the obstacles in proportion to the problem wavelength. In Experiment Three (§6.3), $\omega$ consists of two disjoint triangular scatterers.

In terms of observed error, each value of $\mathbf{d}$ tested gave very similar results, hence we focus here on the case $\mathbf{d}=(1,1)/\sqrt{2}$, which allows some re-reflections between the obstacles and partial illumination of $\Gamma$, see Figures 5LABEL:sub@fig:exp1D, 7LABEL:sub@fig:exp3D and 8LABEL:sub@fig:exp3D.

We now describe the approximation parameters common to all three experiments. To construct the approximation space $V_{N}^{\mathop{\rm HNA}^{*}}(\Gamma,\gamma)$, we first choose $V_{{N_{\Gamma}}}^{\mathop{\rm HNA}}(\Gamma)$ to be the single-mesh approximation space of §4 with $p_{j}=p$ for each side $j=1,\ldots,{\mathcal{N}_{\Gamma}}=3$, reducing the polynomial degree close to the corners of $\Gamma$ in accordance with Remark 4.4, hence $p$ now refers to the polynomial degree on the largest mesh elements. We also remove basis elements close to the corners of the mesh on $\Gamma$ in accordance with Remark 4.2, choosing $\alpha_{j}=\max\{(1+p)/4,2\}$, to improve conditioning of the discrete system (61). A grading parameter of $\sigma=0.15$ is used (as in Hewett et al. (2011), where the rationale for this choice is discussed), with $n_{j}=2p$ layers on each graded mesh, for $j=1,{2,3}$ (hence we may choose the constant from Theorem 4.3 as $c_{j}=2$ ).

Theorem 4.3 ensures that we will observe exponential convergence on $\Gamma$ if the polynomial degree is consistent across the mesh, and Proposition 4.7 ensures that we observe exponential convergence on $\gamma$, if $\gamma$ is analytic. In these numerical experiments we test problems where these two conditions are not met, and encouragingly still observe exponential convergence. As hypothesised by Remark 4.4 and Assumption 4.1, our experiments suggest that our method converges exponentially under conditions much broader than those guaranteed by our theory.

For the standard $hp$-BEM space $V_{{N_{\gamma}}}^{{hp}}(\gamma)$, we use the same parameters $p_{\gamma}=p$, $\sigma$ and $c_{j}$ to grade towards the corners of $\gamma$, so the construction of the mesh on $\gamma$ is much the same as on $\Gamma$. The key difference is that on $\gamma$ every mesh element is sufficiently subdivided to resolve the oscillations. The polynomial degree $p_{j}$ is decreased on smaller elements, as on $\Gamma$, in accordance with Remark 4.2.

Figures 5LABEL:sub@fig:exp1L2, 7LABEL:sub@fig:exp2L2 and 8LABEL:sub@fig:exp3L2 show $L^{2}$ convergence on the boundary $\partial D=\Gamma\cup\gamma$, as $p$ increases, for different values of $k$. The markers correspond to the increasing polynomial degree $p=1,\ldots,7$ and the horizontal axis represents the total number of DOFs $N$, which depends on both $p$ and $k$. The reference solution, denoted $\nu_{N^{*}}$, is computed with $p=8$. Additional checks were performed against a high order standard BEM approximation to validate the reference solution. In each experiment that follows, the increased number of oscillations appears to be handled by the increase in ${N_{\gamma}}$ for each $k$ (here ${N_{\Gamma}}$ remains roughly fixed as $k$ increases, and ${N_{\gamma}}$ increases less than linearly with $k$) with exponential convergence in $p$ observed in each case, as predicted by Theorem 5.3 (for analytic $\gamma$). Given exponential convergence in $L^{2}(\partial D$), corollaries 5.6 and 5.8 are sufficient to guarantee exponential point-wise convergence of the domain approximation (67) and the far-field approximation (68).

Remark 6.1** (Quadrature).**

The integrals in (61) and (62) and the $L^{2}$ norms used to estimate the error in Figures 5LABEL:sub@fig:exp1L2, 7LABEL:sub@fig:exp2L2 and 8LABEL:sub@fig:exp3L2 may be oscillatory and singular. In particular, care must be taken when evaluating the triple integral $\left(\mathcal{A}_{\Gamma\shortrightarrow\gamma}\mathcal{G}_{\gamma\shortrightarrow\Gamma}v,w\right)_{L^{2}(\gamma)}$, which contains a singular oscillatory integrand on elements for which ${\mathcal{G}}_{\gamma\shortrightarrow\Gamma}v$ is supported. Standard composite quadrature routines require a large number of weights and nodes. Hence, at higher frequencies, oscillatory quadrature rules should be used (see Deaño et al. (2018) for a review of such methods), while singular integrals should be computed using a suitable quadrature rule (e.g. Huybrechs & Cools (2009)).

6.1 Experiment one

The configuration tested consists of an equilateral triangle $\Omega$ with perimeter $L_{\Gamma}=6\pi$ and a single small triangular scatterer with perimeter $L_{\gamma}=3\pi/5$, with the obstacles separated by a fixed distance of $\mathop{\rm dist}(\Gamma,\gamma)=\sqrt{3}\pi/5$. The configuration can be seen in 5LABEL:sub@fig:exp1D, which shows the real part of the approximation in the domain (67) for $p=8$.

It follows that there are exactly $k$ wavelengths on each side of $\Gamma$ and $k/10$ on each side of $\gamma$. Experiments were run for $k\in\{20,40,80,160\}$ (so the number of wavelengths across the perimeter ${\partial D}$ ranges from 66 to 528). In Figure 6, we show the real part of the solution $v_{N},(N=1122)$ on $\Gamma$ and $\gamma$, for $k=40$. On $\Gamma$, the first side $(s/(2\pi)\in[0,1])$ is the side in shadow, and the third side $(s/(2\pi)\in[2,3])$ is the illuminated side on the right in Figure 5LABEL:sub@fig:exp1D. On these two sides, the effect of the presence of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} is negligible. However, on the middle side $(s/(2\pi)\in[1,2])$, the effect of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} can clearly be seen.

For a fixed number of DOFs $N$, the $L^{2}(\partial D)$ error is approximately the same for each $k$. For each value of $k$ tested, we achieve approximately $1\%$ relative error with approximately $1000$ DOFs. For $k=160$ the combined boundary $\Gamma\cup\gamma$ is 528 wavelengths long, corresponding to approximately two DOFs per wavelength. This illustrates why the method is particularly well suited to problems with one large polygon (for which the high-frequency asymptotics are well understood), and one (or many) small nearby obstacle(s) on which the high frequency asymptotics do not need to be known.

6.2 Experiment two

Now we test the accuracy of our method as the separation (between the large and small obstacle) shrinks with increasing frequency, keeping all other parameters the same as in Experiment One (§6.1). We choose the separation to be

[TABLE]

as is depicted for $k=20$ in Figure 7LABEL:sub@fig:exp2D. Note the decrease in distance when compared with Figure 5LABEL:sub@fig:exp1D.

Despite the obstacles becoming very close together, with a separation of just $3\pi/160<0.06$ at the highest frequency tested, we observe reassuringly similar $L^{2}(\partial D)$ convergence rates (Figure 5LABEL:sub@fig:exp2L2) to Experiment One (Figure 5LABEL:sub@fig:exp1L2). This should not be unexpected, given that (69) satisfies the separation condition (34). Upon closer inspection, the $L^{2}(\partial D)$ error is actually smaller when the obstacles are closer together, notably the final ($p=7$) data points for $k=40,80,160$.

This experiment demonstrates that our method can be applied to high frequency problems in which the obstacles are very close together. This is particularly encouraging when compared with iterative approaches for multiple scattering, which break down when the obstacles are too close together (as discussed in §1).

6.3 Experiment three

Finally, we apply our method to a problem where the small obstacle consists of two small disjoint triangles $\omega_{\gamma}=\omega_{1}\cup\omega_{2}$. Here we take $\omega_{1}$ to be the smaller triangle from Experiment One (§6.1), translated by $(0,1/2)$, and we take $\omega_{2}$ is the smaller triangle from Experiment One flipped horizontally and translated by $(0,-1/2)$. A key difference when compared with the previous two experiments is that this configuration will induce parabolic trapping. As with the previous experiments we have $\mathop{\rm dist}(\Gamma,\gamma)=\sqrt{3}\pi/5$, although now $L_{\gamma}=6\pi/5$. A consequence of this is that there will be twice as many DOFs in the standard basis $V_{{N_{\gamma}}}^{{hp}}(\gamma)$ than were required for the previous experiments, however with this adjustment we observe similar convergence rates (see Figure 8LABEL:sub@fig:exp3L2). Figure 8LABEL:sub@fig:exp3D shows the configuration and the real part of the domain approximation (67) for $p=8$ and $k=20$.

It is clear from Figure 8LABEL:sub@fig:exp3D that the amplitude reaches four times that of the incoming wave, in the region of trapping between the three triangles. The results of this experiment confirm that our method can be extended to configurations of one large obstacle and multiple small obstacles, and is therefore well-suited for efficient modelling of a wide range of trapping phenomena.

7 Conclusions and further work

For a particular class of multiple scattering configurations, we have presented a numerical method which offers a significant reduction in degrees of freedom required at high frequencies, when compared to standard methods. In particular, our method is most effective when one obstacle is much larger than the others. The theoretical estimates presented in §5 rely on a small number of reasonable assumptions, which we prove to hold under certain conditions. However the numerical results of §6 show exponential convergence and stability with respect to the wavenumber in the broader setting where the small obstacle $\gamma$ is not analytic.

As suggested in Remark 6.1, sophisticated quadrature rules are required in conjunction with the proposed method, but these rules can be difficult to implement for oscillatory and singular double and triple integrals. Alternatively, the approximation space of §4.3 may be implemented as a collocation BEM (following the approach of Gibbs et al. (2019)), which would reduce the dimension of each integral by one, making for easier implementation of oscillatory and singular quadrature rules.

The approach detailed in this paper requires at least one (ideally the largest) of the scatterers to be a convex polygon, but extension of this approach to a far broader class of configurations is possible. The key requirement is that the high frequency asymptotics are understood on {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}, which with further work could instead be, e.g., a two-dimensional screen (Hewett et al., 2015), a non-convex obstacle (Chandler-Wilde et al., 2015), or a penetrable obstacle (Groth et al., 2018). Such extensions would not be trivial, however we believe the framework established in this paper lays appropriate groundwork.

In (Chandler-Wilde et al., 2012a), the HNA method is extended from Dirichlet to impedance boundary conditions. By combining such an extension with the approach taken in the present paper, HNA methods may be designed for multiple obstacles with impedance (or Neumann) boundary conditions.

A final area for future work is the case where $\Gamma\cup\gamma$ is connected, such that $\Gamma$ represents the surface of an obstacle on which an HNA basis can be used, whilst $\gamma$ is the component for which we cannot absorb the high-frequency asymptotics into the approximation space. This extension would require more sophisticated bounds on the operator defined by (23).

Acknowledgements

This work was supported by EPSRC (EP/K000012/1 to S.L., PhD studentship to A.G.). A.G. also acknowledges support from the KU Leuven project C14/15/055. A.M. acknowledges support from GNCS-INDAM and from MIUR through the “Dipartimenti di Eccellenza” Programme (2018-2022) – Dept. of Mathematics, University of Pavia. The authors would also like to thank Prof. Markus Melenk for helpful discussions, and the anonymous referees for insightful comments and suggestions.

Appendix A A coercive multiple scattering formulation

In §5 it was noted that there exists a boundary integral formulation of the BVP (2)–(4) which is coercive (sometimes called $V$-elliptic), provided $|\gamma|$ is of the order of one wavelength. With a coercive formulation, it follows by the Lax–Milgram Theorem that the corresponding discrete problem (equivalent to (61)–(62)) is well posed, on any finite dimensional subspace of $L^{2}(\Gamma\cup\gamma)$. We now present this formulation.

For problems of scattering by a single star-shaped obstacle, it was shown in Spence et al. (2011) that the star combined formulation is coercive for problems on a single star-shaped obstacle. In the thesis Gibbs (2017) this formulation was extended to the constellation combined formulation, where it was shown to be coercive for certain configurations consisting of multiple star-shaped obstacles. We present a version with sharper bounds here, specialising the coercivity result to the case of one large obstacle {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega} and one or many small obstacles {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}. We begin by formally defining the configurations of interest:

Definition A.1** (Star- and constellation-shaped).**

A bounded open set $\Upsilon$ with boundary $\partial\Upsilon$ is star-shaped if there exists $\mathbf{x}^{c}\in\Upsilon$ and a Lipschitz continuous $g:S^{1}\rightarrow\mathbb{R},$ where $S^{1}:=\{\hat{\mathbf{x}}\in\mathbb{R}^{2}:|\hat{\mathbf{x}}|=1\}$, such that $g(\hat{\mathbf{x}})>0$ for all $\hat{\mathbf{x}}\in S^{1}$ with

[TABLE]

Intuitively, this may be interpreted as the following: Given any $\mathbf{x}\in\Upsilon$, one can draw a straight line from $\mathbf{x}^{c}$ to $\mathbf{x}$, without leaving $\Upsilon$.

We say a domain is constellation-shaped if it can be represented as the finite union of multiple star-shaped, pairwise disjoint obstacles. In such a case, for each star-shaped component we denote the above $\mathbf{x}^{c}$ parameter by $\mathbf{x}^{c}_{i}$, where $i$ is the index of that component.

We will use the integral operator

[TABLE]

with the surface gradient operator of the fundamental solution as its kernel:

[TABLE]

where $\Phi_{k}$ is as in (5). Now we define our new BIE:

Definition A.2** (Constellation combined formulation).**

For a constellation-shaped domain $\Upsilon$ with boundary $\partial\Upsilon=\cup_{i=1}^{N_{\Upsilon}}\partial\Upsilon_{i}$, with $\partial\Upsilon_{i}$ the boundary of each star-shaped component, we define the constellation-combined operator $\mathcal{A}_{k}:L^{2}(\partial\Upsilon)\rightarrow L^{2}(\partial\Upsilon)$ as

[TABLE]

where $\mathbf{Z}(\mathbf{x})=\mathbf{x}-\mathbf{x}^{c}_{i}$ (with $\mathbf{x}^{c}_{i}\in\Upsilon_{i}$ chosen as $\mathbf{x}^{c}$ for each star-shaped component in Definition A.1) on $\partial\Upsilon_{i}$, for $i=1,\ldots,N_{\Upsilon}$ and $\hat{\eta}(\mathbf{x}):=k|\mathbf{Z}(\mathbf{x})|+\mathrm{i}/2$. This operator yields an alternative BIE to (7), namely

[TABLE]

where the right-hand side data is

[TABLE]

Invertibility of ${\mathcal{A}}_{k}$ follows by Chandler-Wilde et al. (2012b, Theorem 2.41) and is shown in Gibbs (2017, Theorem 5.6). For single star-shaped obstacles, the following is the key result of Spence et al. (2011).

Theorem A.3**.**

Suppose $\Upsilon$ is star-shaped and ${\mathcal{A}}_{k}$ is defined as in Definition A.2. Then the following coercivity result holds:

[TABLE]

where

[TABLE]

In the thesis Gibbs (2017), the above result was extended to configurations of multiple star-shaped obstacles, under additional geometric constraints. These essentially required the obstacles to be sufficiently far apart, when compared with the wavelength and combined perimeter of the configuration. One way to interpret this is by decomposing ${\mathcal{A}}_{k}$ into block operator form (as in (51)), where each off-diagonal block corresponds to the interaction between two disjoint obstacles, and the diagonal blocks correspond to self interactions. It follows by Theorem A.3 that the diagonal operators will be coercive in a constellation-shaped domain. If the interaction between the obstacles is sufficiently small, then any contribution from the off-diagonal terms will be small, and the full block operator will be coercive. It follows from (71) that the kernel of the integral component

[TABLE]

is

[TABLE]

which simplifies to

[TABLE]

We now consider disjoint, star-shaped boundaries $X$ and $Y$, with $\mathbf{x}\in X$ and $\mathbf{y}\in Y$. We can bound the kernel $K$ by considering Definition A.2 and (71), noting $|\mathbf{Z}(\mathbf{x})|\leq\mathop{\rm diam}(X)$, and upper bounds on the Hankel functions from Chandler-Wilde et al. (2009, (1.22),(1.23))

[TABLE]

It follows by the definition of the operator norm, and the Cauchy–Schwarz inequality (see Gibbs (2017, Lemma 5.13) for a more general result) that for disjoint Lipschitz boundaries $X$ and $Y$,

[TABLE]

The bound (73) quantifies the interaction between two disjoint Lipschitz boundaries $X$ and $Y$. The following theorem exploits this bound, deriving a coercivity result for a subclass of configurations considered in this paper - one large and one (or many) small obstacles.

Theorem A.4**.**

Suppose we have a multiple scattering configuration consisting of one large star-shaped obstacle with boundary $\Gamma$, and ${\mathcal{N}_{\gamma}}$ small star-shaped obstacles $\gamma_{i}$ with boundary $\gamma=\cup_{i}\gamma_{i}$. Suppose further that obstacles are pairwise disjoint, such that the minimum distance between any two obstacles is bounded below by $R>0$. Assuming $|\Gamma|\geq|\gamma|$, if

[TABLE]

then the Constellation Combined operator of Definition A.2 is coercive (i.e. satisfies a bound of the form (72)) with coercivity constant

[TABLE]

Proof A.5**.**

To simplify the notation, we shall write $\|{\mathcal{A}}_{k,Y\shortrightarrow X}\|$ to mean $\|{\mathcal{A}}_{k,Y\shortrightarrow X}\|_{L^{2}(Y)\shortrightarrow L^{2}(X)}$. We begin by decomposing the operator into a sum of operators defined on subsets of $\Gamma\cup\gamma$,

[TABLE]

in which we have split the operator into diagonal and off-diagonal terms

[TABLE]

where we have abused the notation of Definition 4.9, which is used differently here to mean:

[TABLE]

where $\mathbbm{1}_{X}$ is an indicator function, equal to one on $X$ and zero otherwise. The diagonal terms can all be bounded via Theorem A.3, yielding

[TABLE]

We want to find conditions under which the right-hand side of the above inequality is positive, hence we require the negative term to be sufficiently small. We bound these off-diagonal terms

[TABLE]

We now split the above norm on ${\mathcal{A}}_{\text{cross}}$ using the triangle inequality noting the terms in (75), and apply the bound (73) to each component,

[TABLE]

where we have used $|\Gamma|\geq\{|\gamma|,2\mathop{\rm diam}(\Gamma),2\mathop{\rm diam}({\gamma})\}$ to simplify terms. Appealing also to the Cauchy-Schwarz inequality, we can write

[TABLE]

which can be used to simplify (79) to obtain

[TABLE]

Noting (77), we require that

[TABLE]

which is equivalent to the condition (74).

We do not expect the above result to be sharp. A key consequence is the following: if $|\gamma|$ is no more than a fixed fraction of a wavelength, the constellation combined formulation is coercive. We conclude this appendix with bounds on two of the key constants of the Galerkin method as outlined in §5, if the constellation combined formulation is used instead of the standard combined formulation. With the standard formulation, we are unable to bound these constants given current available theory.

Corollary A.6**.**

Suppose we reformulate the Galerkin method of §5 instead using the constellation combined formulation of Definition A.2, and that our scattering configuration \Upsilon={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\cup{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega} satisfies the conditions of Theorem A.4. Then the constants $C_{q}(k)$ and $N_{0}(k)$ of Lemma 5.1, Theorem 5.3 and Corollary 5.6 satisfy

[TABLE]

where $C>0$ is a constant which depends only on the geometry of $\Gamma$ and $\gamma$ and $\alpha_{\Gamma\cup\gamma}$ is the coercivity constant from Theorem A.4.

Proof A.7**.**

Given that the conditions of Lemma 5.1 hold, our formulation is coercive. It follows by the Lax–Milgram Theorem that $N_{0}(k)=1$. It follows by Céa’s Lemma that the quasi-optimality constant is

[TABLE]

The norm in the numerator of (80) is $O(k^{1/2})$ for all $k\geq k_{0}$ (Spence et al., 2011, Theorem 4.2).

Finally, we remark that for a given geometry {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}\cup{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega}, there exists a $k_{1}>0$ such that for all $k\geq k_{1}$, Theorem A.4 cannot guarantee coercivity, and consequentially the statements of Corollary A.6 may not be valid. This is because the negative component of $\alpha_{\Gamma\cup\gamma}$ (as defined in Theorem A.4) will become larger in magnitude as $k$ increases, whilst the positive component remains fixed; we require $\alpha_{\Gamma\cup\gamma}>0$ to ensure coercivity.