A hybrid numerical-asymptotic boundary element method for high frequency scattering by penetrable convex polygons

TL;DR

This paper introduces a hybrid boundary element method combining geometrical optics and diffraction for high frequency scattering by convex polygons, achieving high accuracy with fewer degrees of freedom.

Contribution

The paper presents a novel hybrid numerical-asymptotic boundary element method that efficiently captures high frequency oscillations in scattering problems using a combined geometrical optics and diffraction approach.

Findings

Achieves fixed accuracy with frequency-independent degrees of freedom.

Including diffraction improves accuracy by an order of magnitude.

Effective for acoustic and electromagnetic scattering by convex polygons.

Abstract

We present a novel hybrid numerical-asymptotic boundary element method for high frequency acoustic and electromagnetic scattering by penetrable (dielectric) convex polygons. Our method is based on a standard reformulation of the associated transmission boundary value problem as a direct boundary integral equation for the unknown Cauchy data, but with a nonstandard numerical discretization which efficiently captures the high frequency oscillatory behaviour. The Cauchy data is represented as a sum of the classical geometrical optics approximation, computed by a beam tracing algorithm, plus a contribution due to diffraction, computed by a Galerkin boundary element method using oscillatory basis functions chosen according to the principles of the Geometrical Theory of Diffraction. We demonstrate with a range of numerical experiments that our boundary element method can achieve a fixed accuracy of approximation using only a relatively small, frequency-independent number of degrees of freedom. Moreover, for the scattering scenarios we consider, the inclusion of the diffraction term provides an order of magnitude improvement in accuracy over the geometrical optics approximation alone.

Figures

Tables

Table 1. Table 1: Degrees of freedom per wavelength λ 2 subscript 𝜆 2 \lambda_{2} for the results in Figure 8 . The total number of DOF is 456 for the h ​ p ℎ 𝑝 hp -BEM and 416 for the HNA BEM, independent of k 1 subscript 𝑘 1 k_{1} .

$k_1$	5	10	20	40	80	160
$hp$-BEM $\#\mathrm{DOFper}{\lambda }_2$	10.13	5.07	2.53	1.27	0.63	0.32
HNA $\#\mathrm{DOFper}{\lambda }_2$	9.24	4.62	2.31	1.16	0.58	0.29
$hp$-BEM COND	$6.95\times {10}^1$	$2.54\times {10}^2$	$5.51\times {10}^2$	$1.53\times {10}^3$	$4.94\times {10}^3$	$3.20\times {10}^4$
HNA BEM COND	$1.42\times {10}^8$	$1.55\times {10}^6$	$1.93\times {10}^6$	$4.77\times {10}^5$	$1.46\times {10}^5$	$1.08\times {10}^5$

Table 2. Table 2: Number of DOF in the HNA approximation space for the examples considered in Figure  9 .

${\theta }_i$ \$k_1$	10	20	40	80	160
$\frac{\pi }{2}$ (${𝐝}_1$)	416	416	416	416	416
$\frac{5\pi }{12}$ (${𝐝}_2$)	424	424	424	416	408
$\frac{\pi }{3}$ (${𝐝}_3$)	408	408	408	400	400
$\frac{\pi }{4}$ (${𝐝}_4$)	416	408	408	408	400
$\frac{\pi }{6}$ (${𝐝}_5$)	408	408	408	408	400

Equations

$$k_1>0,k_2>0,\text{ and }\mathrm{Im}\left[k_2\right]\ge 0;$$

$$\mu =k_2/k_1$$

$$\mathrm{Im}\left[\alpha \right]\le 0,\mathrm{Im}\left[\alpha k_2^2\right]\ge 0,$$

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

$$\Delta u_1+k_1^2{u}_1=0,$$

$$\Delta u_2+k_2^2{u}_2=0,$$

$$u_1=u_2\text{mbox}and\frac{\partial u_1}{\partial \mathbf{n}}=\alpha \frac{\partial u_2}{\partial \mathbf{n}},$$

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

$$u_1(\mathbf{x})$$

$$u_2(\mathbf{x})$$

$${\mathcal{S}}_j\varphi (\mathbf{x})={\int }_{\Gamma }{\Phi }_j(\mathbf{x},\mathbf{y})\varphi (\mathbf{y})\mathrm{d}s(\mathbf{y}),{\mathcal{D}}_j\varphi (\mathbf{x})={\int }_{\Gamma }\frac{\partial {\Phi }_j(\mathbf{x},\mathbf{y})}{\partial \mathbf{n}(y)}\varphi (\mathbf{y})\mathrm{d}s(\mathbf{y}).$$

$$u^s(\mathbf{x})\sim \frac{{\mathrm{e}}^{\mathrm{i}\pi /4}}{2\sqrt{2\pi }}\frac{{\mathrm{e}}^{\mathrm{i}{k}_{1}r}}{\sqrt{k_{1}r}}F(\hat{\mathbf{x}}),\text{mbox}asr:=|\mathbf{x}|\to \infty ,$$

$$F(\hat{\mathbf{x}})=-{\int }_{\Gamma }{\mathrm{e}}^{-\mathrm{i}{k}_1\hat{\mathbf{x}}\cdot \mathbf{y}}\left(\mathrm{i}{k}_1(\hat{\mathbf{x}}\cdot \mathbf{n}(\mathbf{y}))u_1(\mathbf{y})+\frac{\partial u_1}{\partial \mathbf{n}}(\mathbf{y})\right)\text{mbox}ds(\mathbf{y}),\hat{\mathbf{x}}\in {\mathbb{S}}^1.$$

$$\mathbf{v}:=\left.\left(\begin{array}[]{c}u_{1}\\ \partial u_{1}/\partial\mathbf{n}\end{array}\right)\right|_{\Gamma}=\left.\left(\begin{array}[]{c}u_{2}\\ \alpha\,\partial u_{2}/\partial\mathbf{n}\end{array}\right)\right|_{\Gamma},$$

$$A\mathbf{v}=\mathbf{f},$$

$$\displaystyle A=\begin{pmatrix}\frac{1+\alpha}{2}I+(\alpha D_{2}-D_{1})&S_{1}-S_{2}\\ \alpha(H_{2}-H_{1})&\frac{1+\alpha}{2}I+(\alpha D_{1}^{\prime}-D_{2}^{\prime})\end{pmatrix},\quad\qquad\mathbf{f}=\left(\begin{array}[]{c}u^{i}\\ \alpha\,\partial u^{i}/\partial\mathbf{n}\end{array}\right),$$

$$S_j\varphi (\mathbf{x}):={\int }_{\Gamma }{\Phi }_j(\mathbf{x},\mathbf{y})\varphi (\mathbf{y})\mathrm{d}s(\mathbf{y}),D_j\varphi (\mathbf{x}):={\int }_{\Gamma }\frac{\partial {\Phi }_j(\mathbf{x},\mathbf{y})}{\partial \mathbf{n}(\mathbf{y})}\varphi (\mathbf{y})\mathrm{d}s(\mathbf{y}),$$

$$D_j^{\prime }\varphi (\mathbf{x}):={\int }_{\Gamma }\frac{\partial {\Phi }_j(\mathbf{x},\mathbf{y})}{\partial \mathbf{n}(\mathbf{x})}\varphi (\mathbf{y})\mathrm{d}s(\mathbf{y}),H_j\varphi (\mathbf{x}):=\frac{\partial }{\partial \mathbf{n}(\mathbf{x})}{\int }_{\Gamma }\frac{\partial {\Phi }_j(\mathbf{x},\mathbf{y})}{\partial \mathbf{n}(\mathbf{y})}\varphi (\mathbf{y})\mathrm{d}s(\mathbf{y}).$$

$$\mathbf{v}(\mathbf{x})\approx {\mathbf{v}}^{GO}(\mathbf{x})+{\mathbf{v}}^d(\mathbf{x}),\mathbf{x}\in \Gamma ,$$

$$\underset{\mathbf{x}}{\max }|a|{\mathrm{e}}^{-E\mathbf{e}\cdot \mathbf{x}}>{\mathrm{tol}}_b,$$

$${\mathbf{d}}^r={\mathbf{d}}^i-2({\mathbf{d}}^i\cdot \mathbf{n})\mathbf{n},{\mathbf{e}}^r={\mathbf{e}}^i-2({\mathbf{e}}^i\cdot \mathbf{n})\mathbf{n},$$

$$D_i^2-E_i^2={\eta }^2-{\xi }^2,D_t^2-E_t^2={\tilde{\eta }}^2-{\tilde{\xi }}^2,D_i{E}_i({\mathbf{d}}^i\cdot {\mathbf{e}}^i)=\eta \xi ,D_t{E}_t({\mathbf{d}}^t\cdot {\mathbf{e}}^t)=\tilde{\eta }\tilde{\xi }.$$

$$D_i{\mathbf{d}}^i\cdot \mathbf{t}=D_t{\mathbf{d}}^t\cdot \mathbf{t},E_i{\mathbf{e}}^i\cdot \mathbf{t}=E_t{\mathbf{e}}^t\cdot \mathbf{t},$$

$$\frac{a^r}{a^i}$$

$$\frac{a^t}{a^i}$$

$$D_t$$

$$E_t$$

$${\mathbf{d}}^t\cdot \mathbf{n}\ge 0,\text{(GO1)},$$

$${\mathbf{e}}^t\cdot \mathbf{n}\ge 0,\text{(GO2)},$$

$$\text{use GO1, unless }|{\mathbf{d}}^t\cdot \mathbf{n}|<{\mathrm{tol}}_{GO},\text{ in which case use GO2}.$$

$$|a|{\mathrm{e}}^{-E\mathbf{e}\cdot {\mathbf{x}}_{\ast }}>{\mathrm{tol}}_{bb},$$

$${\mathbf{v}}^d(\mathbf{x})=\sum _{j=1}^{n_s}{\mathbf{v}}_1^j(\mathbf{x})\exp [\mathrm{i}{k}_1{r}_j(\mathbf{x})]+\sum _{j=1}^{n_s}{\mathbf{v}}_2^j(\mathbf{x})\exp [\mathrm{i}{k}_2{r}_j(\mathbf{x})],\mathbf{x}\in \Gamma ,$$

$$x_0:=0,x_i:={\sigma }^{n-i},i=1,2,\dots ,n,$$

Full text

**Keywords: ** Helmholtz equation, transmission problem, high frequency, boundary element method, Geometrical Theory of Diffraction, acoustic and electromagnetic scattering.

1 Introduction

Scattering of time-harmonic acoustic or electromagnetic waves by penetrable (i.e., dielectric) scatterers arises in numerous physical applications, for example the scattering of light by atmospheric ice crystals, important in determining the earth’s radiation balance in global climate models [3]. When the penetrable scatterer and the exterior domain are both homogeneous, a standard computational approach to solving the associated transmission boundary value problem, at least for low/moderate frequencies, is the boundary element method (BEM) (often called the Method Of Moments in the electromagnetic community), which involves the numerical solution of a system of boundary integral equations (BIEs) on the scatterer boundary [50, 51, 30].

However, conventional BEM implementations based on piecewise polynomial approximation spaces are computationally infeasible in the high frequency regime where the scatterer is large relative to the smallest wavelength in the problem, since they suffer from the well-known limitation (suffered by all conventional numerical approaches including finite element, finite difference and spectral methods) that a fixed number of degrees of freedom is required per wavelength in order to accurately represent the oscillatory solution [44]. At high frequencies the wave field is more naturally described by asymptotic theories such as Geometrical Optics (GO) and the Geometrical Theory of Diffraction (GTD) [38, 8, 37], which represent it in a completely different way, as a superposition of ray fields. In principle, these asymptotic approximations have frequency-independent computational cost. However, they are accurate only at sufficiently high frequencies, and require the solution of certain canonical problems, which are not always solvable in closed form.

In this paper we present a hybrid numerical-asymptotic (HNA) BEM for two-dimensional scattering by penetrable convex polygons. In general terms, the HNA approach aims to develop numerical methods which are computationally feasible over the whole frequency range, by explicitly building features of the high frequency solution behaviour into the numerical approximation space. Specifically, rather than using conventional piecewise polynomial BEM basis functions, one instead uses oscillatory basis functions, whose oscillations match those of the GO and GTD ray fields. The HNA approach (reviewed in [15]) has proven very effective for a number of problems involving impenetrable scatterers, including convex [17, 35, 18], nonconvex [16, 33] and curvilinear [42] polygons, two-dimensional colinear screens [34], smooth convex two-dimensional [22, 13, 25, 24] and three-dimensional [26] scatterers, and three-dimensional rectangular screens [32]. All of the HNA methods mentioned above achieve fixed accuracy with a frequency-independent (or only very modestly growing as the frequency increases) number of degrees of freedom. Moreover, many (in particular [17, 35, 16, 33, 34, 25, 24]) are fully convergent and supported by rigorous numerical analysis.

We believe that the method presented in this paper is the first HNA method for any problem involving a penetrable scatterer. Application of the HNA methodology to penetrable scatterers is nontrivial, primarily because the high frequency asymptotic behaviour is considerably more complicated and less well-understood than in the corresponding impenetrable case. In addition to the phenomena of reflection and diffraction that pertain in the impenetrable case, one also has to take into account the refraction of waves into the interior of the scatterer. This generates a complicated superposition of internally reflected and diffracted wave fields. Furthermore, the exact nature of these wave fields is not fully understood, because many of the relevant GTD canonical scattering problems remain unsolved in closed form. In particular, the key open problem of relevance for scattering by penetrable polygons is diffraction by a penetrable wedge, which has not yet been solved in full generality, despite numerous attempts (see, e.g., [46, 49, 14, 1]). Nonetheless, in spite of these difficulties it was demonstrated in our earlier publication [31] that effective HNA approximation spaces can indeed be developed for penetrable scatterers, by adopting appropriate generalizations of the methods already developed for the impenetrable case.

The two key general ideas explored in [31] are that:

1. (i)

In order to design HNA approximation spaces one needs only *partial *information about the GTD solution (specifically phases, and the location of shadow boundaries), not the full GTD solution. In particular, the lack of analytical expressions for corner diffraction coefficients for penetrable scatterers is not a barrier to the development of HNA methods; 2. (ii)

Even if the high frequency oscillatory behaviour is too complicated to be captured completely by an HNA method (because there are too many phases to consider), incorporating just a small number of oscillatory basis functions capturing the leading order asymptotic behaviour might lead to a significant improvement in performance over conventional purely numerical or purely asymptotic methods.

Specifically, for the case of penetrable convex polygons, the approximation space proposed in [31] represents the unknown Cauchy data in the BIE formulation as the GO field (which can be computed analytically using a beam-tracing method - cf. §3.1 below) plus an HNA component designed to capture the primary, secondary, and all higher-order corner-diffracted fields. This HNA component takes the form of a sum of terms $\mathbf{v}_{1}^{j}{\mathrm{e}}^{{\mathrm{i}}k_{1}r_{j}}$ and $\mathbf{v}_{2}^{j}{\mathrm{e}}^{{\mathrm{i}}k_{2}r_{j}}$ (see equation (25) below), where $k_{1}$ and $k_{2}$ are the exterior and interior wavenumbers, $r_{j}$ is Euclidean distance from the $j$th corner of the polygon, and $\mathbf{v}_{1}^{j}$ and $\mathbf{v}_{2}^{j}$ are slowly-varying amplitudes, which are represented by piecewise polynomials on appropriately graded meshes. A BEM implementation was not presented in [31]; instead the performance of the proposed approximation space was investigated by projecting a highly accurate fully numerical reference solution onto the approximation space and computing the resulting approximation error.

In the current paper we present a Galerkin BEM implementation of a refined version of the HNA approximation space developed in [31], demonstrating that it can be used as the basis of an effective numerical method. Regarding the approximation space itself, compared to [31] we have improved: (i) the robustness of the beam-tracing algorithm used to compute the classical GO approximation (see §3.1.1), leading to improved accuracy for certain incident wave directions; and (ii) the efficiency of our meshing strategy for treating beam boundary discontinuities (see §3.1.2 and §3.2), significantly reducing the number of degrees of freedom required to approximate the diffracted wave fields.

The state of the art for high frequency electromagnetic transmission problems in atmospheric physics appears to be the physical-geometric optics hybrid (PGOH) method of Bi et al. [5, 6]. The PGOH method is an extension to penetrable scatterers of the classical Kirchhoff approximation, in which one takes a boundary integral representation of the solution in the propagation domain (Green’s representation theorem in acoustics; the Stratton-Chu formula in electromagnetics) and replaces the unknown boundary data by its GO approximation (including contributions not just from the incident wave, but also from internally reflected and refracted rays). Our method improves on the PGOH approximation by adding in numerically computed diffracted fields within the integral representation. We shall demonstrate, by a range of numerical experiments that, for all of the examples we consider, our method achieves an order of magnitude improvement in accuracy over the PGOH approach, with only a modest (and frequency-independent) number of degrees of freedom required to compute the additional diffracted contribution.

The outline of the paper is as follows. We begin in §2 by stating the scattering problem and its reformulation as a BIE. In §3 we describe the construction of our HNA approximation space. In §4 we describe the implementation of our approximation space as a Galerkin BEM. We present numerical results demonstrating its efficiency and accuracy for a range of examples with different (convex polygonal) geometries and different physical parameters, including variations in contrast, absorption, and incident angle. In §5 we offer some conclusions and discuss potential future improvements and generalizations of the method.

2 Problem statement

We consider the two-dimensional scattering of a time-harmonic incident plane wave by a penetrable convex polygon. Let $\Omega_{2}\subset\mathbb{R}^{2}$ denote the interior of the polygon (a convex bounded open set), and let $\Omega_{1}:=\mathbb{R}^{2}\backslash\overline{\Omega_{2}}$ denote the exterior unbounded domain. Let $\Gamma=\Gamma_{1}\cup\Gamma_{2}\cup\ldots\Gamma_{n_{s}}$ denote the boundary of the polygon, where $n_{s}$ is the number of sides and $\Gamma_{j}$, $j=1,\ldots,n_{s}$, are the sides of the polygon, labelled in an anti-clockwise direction. The corners of the polygon are similarly labelled $\mathbf{P}_{1},\ldots,\mathbf{P}_{n_{s}}$, with $\Gamma_{j}$, $j=1,\ldots,n_{s}$, being the side between the corners $\mathbf{P}_{j}$ and $\mathbf{P}_{j+1}$ (with the convention that $\mathbf{P}_{n_{s}+1}\equiv\mathbf{P}_{1}$); for an example see Figure 7(a). Let $\mathbf{n}$ denote the outward unit normal to $\Gamma$. Let $k_{1}$ denote the wavenumber in $\Omega_{1}$ and $k_{2}$ the wavenumber in $\Omega_{2}$. We shall assume throughout that

[TABLE]

when ${\rm Im}\left[k_{2}\right]>0$ the scatterer is partially absorbing. For convenience we introduce the notation

[TABLE]

for the (possibly complex) refractive index, and $\lambda_{1}=2\pi/k_{1}$, $\lambda_{2}=2\pi/k_{2}=\lambda_{1}/\mu$ for the exterior and interior wavelengths.

The boundary value problem (BVP) we wish to solve is: given $\alpha\in\mathbb{C}\setminus\{0\}$ satisfying

[TABLE]

and an incident field

[TABLE]

with $\mathbf{d}^{i}\in\mathbb{R}^{2}$ a unit direction vector, determine the total field $u_{1}$ in $\Omega_{1}$ and $u_{2}$ in $\Omega_{2}$ such that

[TABLE]

with the scattered field $u^{s}:=u_{1}-u^{i}$ satisfying the Sommerfeld radiation condition

[TABLE]

The scalar BVP (3)–(7) is a model for both acoustic and electromagnetic scattering problems. In the electromagnetic case, scattering by a (possibly partially absorbing) dielectric polygon is modelled by two BVPs of the form (3)–(7), one for the out-of-plane electric field and one for the out-of-plane magnetic field. The standard transmission boundary conditions for Maxwell’s equations (see, e.g., [7, §1.1]) imply that for the electric field the appropriate choice of $\alpha$ in (6) is $\alpha=1$; for the magnetic field it is $\alpha=1/\mu^{2}=(k_{1}/k_{2})^{2}$. Assuming (1), both of these choices of $\alpha$ satisfy (2).

Under the assumptions (1) and (2) it is well known that the BVP (3)–(7) is uniquely solvable with $u_{1}\in C^{2}(\Omega_{1})\cap H^{1}_{\rm loc}(\Omega_{1})$ and $u_{2}\in C^{2}(\Omega_{2})\cap H^{1}(\Omega_{2})$ (see, e.g., [41, Proposition 2.1 and Corollary 3.4], which follows from results in [20, 52], and also the related result of [45, Corollary 8.5]). We mention also [47], where a particularly simple proof of well-posedness is given for the case $k_{1},k_{2}>0$, along with wavenumber-explicit stability bounds.

Our BEM for solving (3)–(7) is based on a direct BIE formulation. We first observe that if $u_{1}$ and $u_{2}$ satisfy (3)–(7), then a form of Green’s representation theorem holds; i.e., the fields $u_{1}$ and $u_{2}$ in $\Omega_{1}$ and $\Omega_{2}$ can be represented in terms of their Dirichlet and Neumann traces on $\Gamma$. By applying the transmission conditions (6) we can work with the boundary traces of just one of the fields; we choose to work with $u_{1}$, denoting its Dirichlet and Neumann traces on $\Gamma$ by $u_{1}$ and $\partial{u_{1}}/\partial{\mathbf{n}}$ respectively. Then [15, Theorems 2.20 and 2.21]

[TABLE]

where $\mathcal{S}_{j}:H^{-1/2}(\Gamma)\to H^{1}_{\rm loc}(\mathbb{R}^{2}\setminus\Gamma)$ and $\mathcal{D}_{j}:H^{1/2}(\Gamma)\to H^{1}_{\rm loc}(\mathbb{R}^{2}\setminus\Gamma)$, $j=1,2$, are the standard single- and double-layer potentials, defined for sufficiently smooth arguments (for a detailed discussion of the properties of $\mathcal{S}_{j}$ and $\mathcal{D}_{j}$ on Lipschitz domains see [15]) by the integral formulas

[TABLE]

Here $\Phi_{j}(\mathbf{x},\mathbf{y}):=({\mathrm{i}}/4)H_{0}^{(1)}(k_{j}|\mathbf{x}-\mathbf{y}|)$, $j=1,2$, are the fundamental solutions of the Helmholtz equations (4) and (5), respectively, with $H_{0}^{(1)}$ denoting the Hankel function of the first kind of order zero. For later reference we note that the far-field behaviour of the scattered field is

[TABLE]

where $\hat{\mathbf{x}}:=\mathbf{x}/|\mathbf{x}|\in\mathbb{S}^{1}$, the unit circle, and the far-field pattern $F:\mathbb{S}^{1}\to\mathbb{C}$ is given by

[TABLE]

To derive a BIE from (8)–(9) one takes Dirichlet and Neumann traces of both equations and applies the standard jump relations for traces of $\mathcal{S}_{j}$ and $\mathcal{D}_{j}$ (see, e.g., [15, p. 115]). This produces a set of four equations satisfied by the unknown boundary data

[TABLE]

which can be combined in a number of different ways to produce different BIE formulations. The particular formulation we consider involves the BIE

[TABLE]

where

[TABLE]

$I$ is the identity operator and $S_{j}$, $D_{j}$, $D_{j}^{\prime}$, $H_{j}$, $j=1,2$, are, respectively, the single-layer, double-layer, adjoint double-layer, and hypersingular integral operators, defined for sufficiently smooth $\phi$ by

[TABLE]

The operator $A$ in (14) is bounded and invertible as a mapping from $L^{2}(\Gamma)\times L^{2}(\Gamma)\to L^{2}(\Gamma)\times L^{2}(\Gamma)$; this can be proved by a straightforward modification of the arguments presented in [52, Proof of Theorem 7.2] (for details see [29, §2.7]). That the solution of the BVP is sufficiently smooth to work in this setting follows from the boundedness and invertibility of $A$ as a mapping from $H^{1}(\Gamma)\times L^{2}(\Gamma)\to H^{1}(\Gamma)\times L^{2}(\Gamma)$ (again, see [52, 29]) and the fact that $u^{i}$ is smooth, so that $\mathbf{f}\in H^{1}(\Gamma)\times L^{2}(\Gamma)$. In fact, by standard elliptic regularity results (see, e.g., [28]) the solution $\mathbf{v}$ is infinitely smooth on each of the sides of the polygon. Singular behaviour of $\mathbf{v}$ at the corners of the polygon will be captured in our BEM using mesh refinement.

Our choice of the particular formulation (11), which is similar to that used in [19, §3.8] (where only smooth scatterers are considered) and also to that used in [52] (albeit for an indirect method, in which the unknowns are non-physical “densities”, rather than the boundary data itself), was made because the cancellation of the strong singularities between the two hypersingular operators in the term $H_{2}-H_{1}$ makes implementation particularly simple. However, we emphasise that our HNA approximation space (described in the next section) can be applied in the context of any direct BIE formulation in which the BIE solution is some combination of the physical unknowns $u_{1}$ and $\partial{u_{1}}/\partial{\mathbf{n}}$, see, e.g., [20, 40, 41, 36, 23].

We end this section with a further comment on the choice of BIE formulation for high frequency transmission problems. When using conventional (i.e., non-HNA) approximation spaces, it is well-known that conditioning issues can severely limit the performance of iterative methods that are necessary for solving the large linear systems that arise. Because of this, the development of new well-conditioned BIE formulations is a highly active area of current research; in particular we mention the indirect formulations recently proposed and analysed in [10, 11, 12], which offer favourable spectral properties compared to other existing formulations. (The recent paper [23] describes a direct counterpart with similar properties.) By contrast, we emphasise that the systems arising from our HNA approach are small enough that they can be solved directly rather than iteratively. Indeed, the whole aim of HNA methods is to obtain small linear systems whose size does not grow with respect to frequency, hence reducing the importance of conditioning.

3 HNA approximation space

Our BEM solves the integral equation (11) by the Galerkin method, using an HNA approximation space specially designed to capture the high frequency behaviour of the solution $\mathbf{v}=(u_{1},\partial{u_{1}}/\partial{\mathbf{n}})^{T}$. Specifically, our method is based on the ansatz

[TABLE]

where $\mathbf{v}^{GO}=(u^{GO},\partial{u^{GO}}/\partial{\mathbf{n}})$ is the GO approximation (computed using the beam tracing algorithm described in §3.1), and $\mathbf{v}^{d}=(u^{d},\partial{u^{d}}/\partial{\mathbf{n}})$ is a contribution from diffracted fields (computed using the HNA approximation strategy described in §3.2). Before spelling out in more detail the practicalities of how $\mathbf{v}^{GO}$ and $\mathbf{v}^{d}$ are computed, we pause to briefly explain the origin of (15) and discuss the accuracy of approximation one might expect it to provide.

As was mentioned in the introduction, the high frequency asymptotic theory for transmission problems is not fully understood, because the details of the GTD for penetrable scatterers, in particular the solution of the canonical penetrable wedge problem, are yet to be fully worked out. However, by analogy with the impenetrable case, one expects the leading order high frequency behaviour of $\mathbf{v}$ to be given by the GO approximation $\mathbf{v}^{GO}$, which consists of the incident wave, along with reflected and transmitted waves satisfying Snell’s law and the Fresnel formulas. For penetrable polygons irradiated by a plane wave, neighbouring parallel rays incident on the same side of the polygon remain parallel under reflection/transmission, and hence $\mathbf{v}^{GO}$ can be computed efficiently using a beam tracing algorithm, as we describe in §3.1. As mentioned in §1, the field approximation obtained by substituting this GO approximation into the integral representation formulas (8) and (9) is referred to as the PGOH method in the atmospheric physics literature, and can be interpreted as a generalization of the classical Kirchhoff approximation for impenetrable scatterers.

As the next-order correction to the GO approximation, we expect diffracted wave fields emanating radially outwards from each of the corners of the polygon into the exterior domain $\Omega_{1}$ (with wavenumber $k_{1}$) and inwards into the interior domain $\Omega_{2}$ (with wavenumber $k_{2}$). It is these diffracted fields, denoted $\mathbf{v}^{d}$, that we compute numerically using our BEM, using the HNA approximation strategy described in §3.2; see in particular the ansatz in equation (25).

The contribution made by the diffracted term $\mathbf{v}^{d}$ is illustrated in Figures 1 and 2, which relate to scattering by an equilateral triangle of refractive index $\mu=1.5+0.003125{\mathrm{i}}$ (for full details of the setup see §4). Figure 1(a) shows the real part of the total field computed using our HNA BEM, i.e., first calculating the approximation $\mathbf{v}\approx\mathbf{v}^{GO}+\mathbf{v}^{d}$ and then substituting this into the representation formulas (8) and (9). Figure 1(b) shows the contribution made by the diffracted term $\mathbf{v}^{d}$ alone; here one can clearly see the circular diffracted waves emanating from the vertices, and the complicated interference pattern they produce inside the scatterer.

Figure 2 contains plots of the solution $u_{1}$ on $\Gamma$, i.e., the first component of the vector $\mathbf{v}$. The upper panels show the real part of $u_{1}$ on $\Gamma$, computed both using a highly accurate reference solution (described in §4) and using the GO approximation $u^{GO}$. The lower panels show the difference $u_{1}-u^{GO}$; this remainder is what we approximate using our diffracted term $u^{d}$. A comparison of the left-hand panels ($k_{1}=10$) with the right-hand panels ($k_{1}=40$), suggests that as the wavenumber $k_{1}$ increases the error $u_{1}-u^{GO}$ decreases in magnitude, at least away from the corners of the polygon. (Similar behaviour is observed for the GO approximation to $\partial{u_{1}}/\partial{\mathbf{n}}$, not shown here.) This is in line with what one should expect from the high frequency GO approximation, i.e., that it becomes more accurate with increasing $k_{1}$, but that it does not capture any effects due to the corner diffraction.

We also expect the high frequency solution behaviour to feature other higher-order effects, such as so-called *lateral *waves (sometimes known as head or bow waves in related contexts), and the internal re-reflections of the diffracted and lateral wave fields. (For a more detailed discussion see [31, §3.2.1], and for a schematic illustration of some of these phenomena see Figure 3.) However, to limit the implementational complexity of our method these higher-order effects are not incorporated into our approximation space. Hence, in contrast to the HNA approximation spaces developed in [17, 18, 35] for impenetrable convex polygons, our HNA approximation space does not capture all of the oscillatory solution behaviour. Accordingly, our method remains “asymptotic” in nature, with an inherent frequency-dependent error which decreases as the frequency tends to infinity. 333Many HNA methods share this same “asymptotic” nature. In particular we note that the methods of [22, 13, 26] for smooth convex impenetrable scatterers all incur an exponentially small frequency-dependent error, since they approximate by zero in the deep shadow region, neglecting the exponentially small creeping wave fields. In contrast, the method of [27] and the analysis of [2] do explicitly incorporate creeping wave fields, and the methods of [25, 24], while not fully capturing creeping wave fields, include enough degrees of freedom in the deep shadow to be fully convergent, without incurring a systematic frequency-dependent error. Nonetheless, our numerical results in §4 show that our method achieves a significant improvement over using GO alone. As we shall see in §4, our neglect of the higher-order terms (lateral waves, and the re-reflections of the diffracted and lateral waves) has less impact when the scatterer is partially absorbing (${\rm Im}\left[k_{2}\right]>0$), since then the neglected terms are attenuated and play less of a role.

We now describe the computation of each of the two terms $\mathbf{v}^{GO}$ and $\mathbf{v}^{d}$ in more detail.

3.1 Computation of $\mathbf{v}^{GO}$

Our GO approximation $\mathbf{v}^{GO}$ in (15) is computed using a beam-tracing algorithm (BTA), similar to those presented in [9] and [6] for the analogous 3D problem. For full details of our BTA we refer the reader to [31, §3.1]; here we simply sketch the basic algorithm and point out some modifications in our current implementation compared to that described in [31].

According to the principles of GO, when the incident plane wave impinges on a side of the polygon $\Omega_{2}$, it generates a beam of reflected rays propagating back into the exterior domain $\Omega_{1}$, and a beam of transmitted rays propagating into the interior domain $\Omega_{2}$. The transmitted beam then generates a sequence of higher-order internally reflected beams and externally transmitted beams, as illustrated in Figure 4. Each beam is described by

1. (i)

a plane wave $a{\mathrm{e}}^{{\mathrm{i}}(D\mathbf{d}+{\mathrm{i}}E\mathbf{e})\cdot\mathbf{x}}$, with amplitude $a$, unit propagation and decay direction vectors $\mathbf{d}$ and $\mathbf{e}$, and (wavenumber-dependent) propagation/decay constants $D>0$ and $E\geq 0$, all determined by the well-known laws of reflection and refraction at an interface between two homogeneous media (Snell’s law (18) and the Fresnel formulas (19)–(20) below), and 2. (ii)

(up to) two beam boundaries outside of which we cut off the plane wave beam sharply to zero.

In principle there is an infinite sequence of internal reflections to consider, but for a practical implementation we need to truncate the sequence somehow. Our truncation criterion is as follows. Suppose that a beam $a{\mathrm{e}}^{{\mathrm{i}}(D\mathbf{d}+{\mathrm{i}}E\mathbf{e})\cdot\mathbf{x}}$ (with associated beam boundaries) illuminates a portion of a side $\Gamma_{j}$. Then we include its contribution to $\mathbf{v}^{GO}$ on $\Gamma_{j}$ if and only if

[TABLE]

where the maximum is taken over all $\mathbf{x}\in\Gamma_{j}$ which lie between the two beam boundaries of the beam in question, and ${\rm tol}_{b}$ is a user-specified tolerance. If a beam fails this test, we exclude it from the calculation of $\mathbf{v}^{GO}$ on $\Gamma_{j}$ and do not track any higher-order beams generated by the reflection of that beam from the side $\Gamma_{j}$. The exact number of beams traced depends upon ${\rm tol}_{b}$, frequency, absorption, geometry, and incident angle. Generally speaking, higher absorption typically leads to fewer beams, whereas a more complex geometry (more sides) leads to more beams. For the numerical results in §4 we used ${\rm tol}_{b}=0.005$; this was found to guarantee convergence of the BTA to within a relative accuracy of around $0.1\%$ for all the examples considered, using between $5$ and $30$ beams, depending on the example. We emphasize that tracing each beam requires us to trace just two geometric rays - the two rays forming the beam boundaries. As a result, in all of the examples we consider, the computational cost of the BTA (both in terms of memory and computation time) is negligible within the overall cost of the HNA BEM, even when a relatively large number of beams is required (such as for the hexagon in §4.8).

3.1.1 Sign choice in the Fresnel formulas

Regarding point (i) above, we remark that the laws of reflection and refraction are completely classical in the case when both media are non-absorbing (see, e.g., [7]). However, when one of the media is absorbing (here, this corresponds to ${\rm Im}\left[k_{2}\right]>0$), there is some debate in the literature about how to make a particular sign choice in the Fresnel formulas. We reviewed this issue in [31, Appendix A], providing an analysis of the general interface problem (with up to two absorbing media) as well as numerous references to relevant literature. But since the publication of [31], further investigations have led us to refine our sign choice rule, as we now explain.

Let us consider the interface problem between a medium with wavenumber $k=\eta+{\mathrm{i}}\xi$ and a second medium with wavenumber $\tilde{k}=\tilde{\eta}+{\mathrm{i}}\tilde{\xi}$, where $\eta,\tilde{\eta}>0$, $\xi,\tilde{\xi}\geq 0$. Let $\mathbf{t}$ and $\mathbf{n}$ be unit tangent and normal vectors on the interface (assumed to contain the origin), with $\mathbf{n}$ pointing into the second medium. An incident plane wave $u^{i}=a^{i}\exp[{\mathrm{i}}(D_{i}\mathbf{d}^{i}+{\mathrm{i}}E_{i}\mathbf{e}^{i})\cdot\mathbf{x}]$ propagating in the first medium generates a reflected wave $u^{r}=a^{r}\exp[{\mathrm{i}}(D_{i}\mathbf{d}^{r}+{\mathrm{i}}E_{i}\mathbf{e}^{r})\cdot\mathbf{x}]$ (also propagating in the first medium), and a transmitted wave $u^{t}=a^{t}\exp[{\mathrm{i}}(D_{t}\mathbf{d}^{t}+{\mathrm{i}}E_{t}\mathbf{e}^{t})\cdot\mathbf{x}]$ (propagating in the second medium). According to the law of reflection the vectors $\mathbf{d}^{r}$ and $\mathbf{e}^{r}$ should satisfy

[TABLE]

and for $u^{i}$, $u^{r}$ and $u^{t}$ to be solutions of the relevant Helmholtz equations we require

[TABLE]

Assuming the boundary conditions $u^{i}+u^{r}=u^{t}$ and $\partial{u^{i}}/\partial{\mathbf{n}}+\partial{u^{r}}/\partial{\mathbf{n}}=\alpha\partial{u^{t}}/\partial{\mathbf{n}}$ on the interface, for some $\alpha\in\mathbb{C}$, one can derive Snell’s law

[TABLE]

and the Fresnel formulas

[TABLE]

in which

[TABLE]

with ${\tilde{D}_{i}}:=D_{i}\mathbf{d}^{i}\cdot\mathbf{t}$ and ${\tilde{E}_{i}}:=E_{i}\mathbf{e}^{i}\cdot\mathbf{t}$ (see the appendix of [31] for more detail). Snell’s law and the Fresnel formulas respectively specify the tangential components of the transmitted vectors $\mathbf{d}^{t}$ and $\mathbf{e}^{t}$, and the amplitudes $a^{r}$ and $a^{t}$. Since $\mathbf{d}^{t}$ and $\mathbf{e}^{t}$ have unit length, their normal components are thus specified up to sign. However, to determine these signs and close the problem one must impose an additional constraint based on some physical argument. Two common constraints require either that

[TABLE]

which prohibits energy flow from the second medium to the first, or that

[TABLE]

which prohibits exponential growth in the second medium.

For some configurations these two choices are compatible, and can both be satisfied simultaneously. As an example we consider the phenomenon of total internal reflection, which occurs when both media are non-absorbing ($\xi=\tilde{\xi}=0$), the second medium has a faster propagation speed than the first (i.e., $0<\tilde{\eta}<\eta$) and the incident direction is close to grazing (i.e., $|\mathbf{d}^{i}\cdot\mathbf{t}|$ is sufficiently small). In this case $\mathbf{d}^{t}\cdot\mathbf{n}=0$, so GO1 is automatically satisfied. But $\mathbf{e}^{t}$ can be either $\mathbf{n}$ or $-\mathbf{n}$. Enforcing GO2 selects the first option, and we obtain the configuration shown in Figure 5(a).

However, when one or both of the media are absorbing, there exist configurations for which GO1 and GO2 are incompatible. (Note that $\mathbf{d}^{t}$ and $\mathbf{e}^{t}$ are coupled through the final equation in (17) so the signs of $\mathbf{d}^{t}\cdot\mathbf{n}$ and $\mathbf{e}^{t}\cdot\mathbf{n}$ cannot in general be chosen independently.) For the case when just one of the media is absorbing (as in the current paper), the only situation in which this incompatibility arises is when the first medium is absorbing ($\xi>0$) and the second is non-absorbing ($\tilde{\xi}=0$), and the incident vectors $\mathbf{d}^{i}$ and $\mathbf{e}^{i}$ both point in the same direction relative to $\mathbf{t}$ (i.e., $(\mathbf{d}^{i}\cdot\mathbf{t})(\mathbf{e}^{i}\cdot\mathbf{t})>0$). In this case Snell’s law implies that $(\mathbf{d}^{t}\cdot\mathbf{n})(\mathbf{e}^{t}\cdot\mathbf{n})<0$, so that (21) and (22) cannot both hold simultaneously. The two possibilities are illustrated in Figure 5(b).

In the context of our problem, this situation can arise when the scatterer is absorbing (${\rm Im}\left[k_{2}\right]>0$) and a beam propagating in $\Omega_{2}$ is incident on $\Gamma$, generating an internally reflected beam and an externally transmitted beam. We expect that the correct choice of GO1 or GO2 in this context could in principle be resolved by matching the local asymptotic behaviour near the interface with the global behaviour of the scattered field. However, due to the lack of a convenient expression for the field near the corners of the polygon we have been unable to carry out such an analysis. Instead we have developed a heuristic approach for determining which choice to make, based on the results of numerical experiments. In these experiments (detailed more fully in [29, §4.3.4]) the GO approximation $\mathbf{v}^{GO}$, computed using first GO1 and then GO2, was compared to a conventional BEM reference solution, for a particular geometrical configuration chosen specifically to isolate the behaviour in question. We found that GO1 gave a more accurate approximation than GO2, except when $\mathbf{d}^{t}$ was almost parallel to the interface. On the basis of these results we adopt the following strategy: whenever GO1 and GO2 are incompatible (as described above),

[TABLE]

Here ${\rm tol}_{GO}$ is a user-defined tolerance; for the numerical results in §4 we used the value ${\rm tol}_{GO}=0.01$. (In our previous paper [31] we used GO1 throughout, which corresponds to setting ${\rm tol}_{GO}=0$. The significant improvement in accuracy resulting from the apparently minor change to ${\rm tol}_{GO}=0.01$, enabling GO2 to be used where appropriate, is described in [29, §4.3.4].)

3.1.2 Beam boundary discontinuities

Regarding point (ii) above, we note that cutting off each beam sharply across its two beam boundaries introduces artificial discontinuities in our approximation (15) to the function $\mathbf{v}$ at the “beam boundary points” where the beam boundaries intersect $\Gamma$ (see Figure 4). If these discontinuities are sufficiently large in magnitude they can cause difficulties for our numerical calculation of the diffracted component $\mathbf{v}^{d}$. We discuss this issue in more detail in the next section, and for now we simply introduce some terminology. Suppose that a beam $a{\mathrm{e}}^{{\mathrm{i}}(D\mathbf{d}+{\mathrm{i}}E\mathbf{e})\cdot\mathbf{x}}$ (with associated beam boundaries) impinges on a side $\Gamma_{j}$, with one of its beam boundaries intersecting $\Gamma_{j}$ at the point $\mathbf{x}_{*}\in\Gamma_{j}$. Then we designate the resulting beam boundary discontinuity as “strong” if

[TABLE]

where ${\rm tol}_{bb}$ is a user-specified tolerance.

3.2 Computation of $\mathbf{v}^{d}$

Our diffraction contribution $\mathbf{v}^{d}$ in (15) takes the form

[TABLE]

where, for $j=1,\ldots,n_{s}$, $r_{j}(\mathbf{x})=|\mathbf{x}-\mathbf{P}_{j}|$ is the Euclidean distance between $\mathbf{x}\in\Gamma$ and the corner $\mathbf{P}_{j}$, and the amplitudes $\mathbf{v}_{l}^{j}$, $l=1,2$, are computed numerically using our Galerkin BEM (see (30) below), as (vector-valued) piecewise polynomials on appropriate overlapping meshes. (Recall that since $\mathbf{v}=(u_{1},\partial{u_{1}}/\partial{\mathbf{n}})$ is vector-valued, so are the amplitudes $\mathbf{v}_{l}^{j}$.)

A typical term $\mathbf{v}_{1}^{j}(\mathbf{x})\exp[{\mathrm{i}}k_{1}r_{j}(\mathbf{x})]$ in the first sum in (25) corresponds to a diffracted wave emanating from the corner $\mathbf{P}_{j}$ and propagating with wavenumber $k_{1}$ in the exterior domain $\Omega_{1}$. The convexity of the polygon $\Omega_{2}$ implies that this wave only impinges on the sides adjacent to the corner $\mathbf{P}_{j}$, namely $\Gamma_{j}$ and $\Gamma_{j-1}$ (with the convention that $\Gamma_{0}\equiv\Gamma_{n_{s}}$). Hence we force $\mathbf{v}_{1}^{j}(\mathbf{x})=\boldsymbol{0}$ for $\mathbf{x}\in\Gamma\setminus(\Gamma_{j}\cup\Gamma_{j-1})$. On each of $\Gamma_{j}$ and $\Gamma_{j-1}$, to capture the singular behaviour at the corner $\mathbf{P}_{j}$ we take $\mathbf{v}_{1}^{j}$ to be a piecewise polynomial on a mesh graded geometrically towards $\mathbf{P}_{j}$ (more details of these meshes are given below).

A typical term $\mathbf{v}_{2}^{j}(\mathbf{x})\exp[{\mathrm{i}}k_{2}r_{j}(\mathbf{x})]$ in the second sum in (25) corresponds to a diffracted wave emanating from the corner $\mathbf{P}_{j}$ and propagating with wavenumber $k_{2}$ in the interior domain $\Omega_{2}$. By the convexity of $\Omega_{2}$, this wave impinges on all sides of the polygon. Hence $\mathbf{v}_{2}^{j}(\mathbf{x})$ is supported on the whole of $\Gamma$. As for $\mathbf{v}_{1}^{j}$, on each of $\Gamma_{j}$ and $\Gamma_{j-1}$ we take $\mathbf{v}_{2}^{j}$ to be a piecewise polynomial on a mesh graded geometrically towards $\mathbf{P}_{j}$ (again, more details of these meshes are given below).

On $\Gamma\setminus(\Gamma_{j}\cup\Gamma_{j-1})$ we take $\mathbf{v}_{2}^{j}$ to be a piecewise polynomial of maximum degree $p\in\mathbb{N}_{0}$ on a mesh of $\Gamma\setminus(\Gamma_{j}\cup\Gamma_{j-1})$ constructed in the following way. We start with a mesh of $n_{s}-2$ elements, consisting of the $n_{s}-2$ sides $\Gamma_{l}$, $l\not\in\{j,j-1\}$ (with $n_{s}-1$ mesh points at the corners $\mathbf{P}_{l}$, $l\neq j$). We then insert mesh points at any beam boundary points $\mathbf{x}\in\Gamma\setminus(\Gamma_{j}\cup\Gamma_{j-1})$ associated with “strong” beam boundaries (as defined in (24)) in $\mathbf{v}^{GO}$, provided that the beam boundaries in question emanated from the corner $\mathbf{P}_{j}$. (Here we note another difference between the current algorithm and that described in [31], where mesh points were introduced at all strong beam boundaries, not just those emanating from corner $\mathbf{P}_{j}$, leading to redundancy in the approximation space due to the same mesh points being added to multiple overlapping meshes. The current algorithm avoids this redundancy, giving a significant reduction in the total number of fewer degrees of freedom, without compromising on accuracy.)

To clarify the mesh refinement procedure, we consider the contribution of the beams illustrated in Figure 4. Let us assume that all of the beam boundary discontinuities in Figure 4 are strong (i.e., satisfy (24)). Then for the primary transmitted beam emanating from side $\Gamma_{2}$ in Figure 4(a), we would insert a mesh point at $\mathbf{x}_{5}$ in the mesh for $\mathbf{v}_{2}^{2}$ (associated with diffraction from corner $\mathbf{P}_{2}$), and at $\mathbf{x}_{4}$ in the mesh for $\mathbf{v}_{2}^{3}$ (associated with diffraction from corner $\mathbf{P}_{3}$). For the internally reflected beams emanating from side $\Gamma_{4}$ and $\Gamma_{5}$ in Figure 4(b), we would insert a mesh point at $\mathbf{x}_{3}$ in the mesh for $\mathbf{v}_{2}^{5}$ (associated with diffraction from corner $\mathbf{P}_{5}$). But we would not insert mesh points at $\mathbf{x}_{2}$ and $\mathbf{x}_{6}$ in the mesh for any of the amplitudes $\mathbf{v}_{2}^{j}$, $j=1,\ldots,n_{s}$, as the points $\mathbf{x}_{2}$ and $\mathbf{x}_{6}$ correspond to beam boundaries that are not directly associated with any corner $\mathbf{P}_{j}$. (Rather, they arise by the internal reflection of the two beam boundaries in Figure 4(a), and are associated with the higher-order diffracted-reflected fields, which are not captured in our approximation space.)

To describe in more detail the geometrically graded meshes we use, we first consider a mesh on the interval $[0,1]$, refined towards [math]. Given $n\geq 1$ (the number of layers in the mesh) we let $G_{n}$ denote the set of meshpoints $\{x_{i}\}_{i=0}^{n}$ defined by

[TABLE]

where $0<\sigma<1$ is a grading parameter (we discuss the choice of $\sigma$ in §4.1). We then define a space of piecewise polynomials on the mesh $G_{n}$, for a given degree vector $\mathbf{p}\in(\mathbb{N}_{0})^{n}$, by

[TABLE]

For reasons of efficiency and conditioning it is common to decrease the degree of the approximating polynomials towards the point of refinement. Specifically, given a maximum polynomial degree $p\in\mathbb{N}_{0}$, we use a “linear slope” degree vector $\mathbf{p}$ with

[TABLE]

We adopt an “$hp$” approximation strategy, in which the number of layers $n$ in each geometric mesh grows as the maximum polynomial degree $p$ is increased. Specifically, we take

[TABLE]

where $C_{np}>0$ is a user-specified constant. The geometric meshes in our approximation space are constructed from this basic building block by straightforward coordinate transformations. Explicitly, for each $j=1,\ldots,n_{s}$ the amplitudes $\mathbf{v}_{1}^{j}$ and $\mathbf{v}_{2}^{j}$ are approximated on the side $\Gamma_{j}$ using the mesh $\mathbf{P}_{j}+G_{n}*(\mathbf{P}_{j+1}-\mathbf{P}_{j})$, and on the side $\Gamma_{j-1}$ using the mesh $\mathbf{P}_{j}+G_{n}*(\mathbf{P}_{j-1}-\mathbf{P}_{j})$.

The above construction constrains $\mathbf{v}^{d}$ to lie in a particular finite-dimensional approximation space $V\subset L^{2}(\Gamma)\times L^{2}(\Gamma)$ whose dimension $N$ (the total number of degrees of freedom) is given by

[TABLE]

where $n_{bb}$ is the total number of strong beam boundary points at which we insert extra mesh points. Note that, by (28), $N$ grows approximately in proportion to $p^{2}$ as $p$ increases.

3.3 Numerical best approximation error analysis

The HNA approximation strategy described above has been extensively tested in [29] and [31] using a numerical best approximation error analysis. For a number of scattering scenarios, highly accurate reference solutions, denoted here by $\mathbf{v}_{\rm ref}$, were computed using a conventional BEM (see §4 for more details). The corresponding GO approximations $\mathbf{v}^{GO}$ were calculated using the BTA, and numerical best approximations to $\mathbf{v}_{\rm ref}-\mathbf{v}^{GO}$ (which we identify with $\mathbf{v}^{d}$, see (15)) from the HNA approximation space $V$ were then computed using orthogonal projection in $L^{2}(\Gamma)\times L^{2}(\Gamma)$ (via the least squares approach detailed in [31]), and the corresponding errors recorded.

A typical sample of the results obtained is presented in Figure 6, which plots the numerical best approximation error against increasing polynomial degree. For reference we also show the error for GO alone, as the left-most data point. All of the curves in Figure 6 exhibit the same qualitative behaviour. At first the error decays exponentially, in an essentially frequency-independent way; in this initial phase, increasing $p$ leads to an improvement in the accuracy of our HNA approximation $\mathbf{v}^{d}$ to the diffracted wave fields. Then, when the accuracy of $\mathbf{v}^{d}$ reaches the level of the asymptotic error associated with our neglect of the higher-order asymptotic phenomena (the lateral waves and the internal re-reflections of the diffracted and lateral waves), the error levels out. Increasing $p$ no longer leads to a significant improvement in accuracy, because our HNA approximation space for $\mathbf{v}^{d}$ does not contain the correct oscillatory basis functions needed to capture the neglected waves. (Exponential convergence will in theory resume once the number of degrees of freedom is sufficiently large to resolve the missing oscillations, but this is not the regime in which HNA methods are designed to operate.) We note that the magnitude of the neglected waves, and hence the level at which the error stagnates, decreases as $k_{1}$ increases, and also as the absorption ${\rm Im}\left[\mu\right]$ increases, in accordance with the asymptotic theory (see the discussion before §3.1).

4 Galerkin BEM and numerical results

In §3 we proposed an ansatz $\mathbf{v}\approx\mathbf{v}^{GO}+\mathbf{v}^{d}$ (see (15)) for the solution $\mathbf{v}=(u_{1},\partial{u_{1}}/\partial{\mathbf{n}})^{T}\in L^{2}(\Gamma)\times L^{2}(\Gamma)$ of the integral equation (11), describing a beam tracing algorithm for computing the GO term $\mathbf{v}^{GO}$ and an HNA approximation space $V\subset L^{2}(\Gamma)\times L^{2}(\Gamma)$ for the diffraction term $\mathbf{v}^{d}$. In this section we describe our BEM implementation and present a range of numerical results demonstrating its performance. Our BEM selects a particular element $\mathbf{v}^{d}\in V$ by applying the Galerkin method to the integral equation (11), rewritten using the decomposition (15), with $\mathbf{v}^{d}$ as the unknown. That is, we compute $\mathbf{v}^{d}\in V$ satisfying

[TABLE]

where $\langle\cdot,\cdot\rangle$ is the usual inner product on $L^{2}(\Gamma)\times L^{2}(\Gamma)$.

4.1 Implementation details

Choosing a basis for $V$ reduces (30) to the solution of a linear system for the basis coefficients of $\mathbf{v}^{d}$. In our implementation we represent each of the two components of $\mathbf{v}^{d}$ using scalar basis functions of the form $L_{m}{\mathrm{e}}^{{\mathrm{i}}k_{l}r_{j}}$ (see (25)), where $l\in\{1,2\}$, $r_{j}$ is the distance from vertex $\mathbf{P}_{j}$, for some $j\in\{1,\ldots,n_{s}\}$, and $L_{m}$ is a Legendre polynomial of degree $m\in\{0,\ldots,p\}$. (Note that each basis function is supported on some subinterval of one of the sides of the polygon.) To improve conditioning of the resulting linear system we found it beneficial to scale each basis function so that its $L^{2}(\Gamma)$ norm is equal to one.

By design, the linear systems produced by the HNA BEM are small compared to those arising from conventional methods and hence can be solved using direct methods, without the need for iterative solvers. Assembly of the linear system requires the evaluation of integrals which are potentially both oscillatory and singular. We evaluate these integrals to high precision (better than $10^{-8}$ relative accuracy) with composite Gauss quadrature using geometric grading and the classical Duffy transformation - for details see [29]. For a faster implementation one could employ oscillatory quadrature techniques, as discussed in [15, 21, 39, 34, 35]. But our focus in the current paper is on demonstrating the efficiency of the HNA approximation space in terms of the number of degrees of freedom, rather than on fast implementations, and so we postpone further discussion of this important issue to future work.

In all of our experiments we use the parameter choices:

[TABLE]

and the fixed maximum polynomial degree

[TABLE]

In accordance with the discussion in §3.3, choosing a larger value of $p$ than $3$ in general leads to improved accuracy, up to the point at which the underlying asymptotic error is reached. However, through extensive numerical testing (detailed in [29]) the choices (31)–(35) were found to give a reasonable compromise between computational effort and solution accuracy for the range of problems considered here.

For the construction of our geometric meshes (recall (26)) we use grading parameters $\sigma_{1}=0.17$ for the $\mathbf{v}_{1}^{j}$ meshes and $\sigma_{2}=0.15$ for the $\mathbf{v}_{2}^{j}$ meshes. Using the same grading parameter for both the $\mathbf{v}_{1}^{j}$ and $\mathbf{v}_{2}^{j}$ meshes was found to cause ill-conditioning when $k_{1}$ is small, or when the meshes are particularly heavily refined, because when the diameter of the smallest elements in the meshes is comparable to the wavelength(s), the polynomial approximants are capable of resolving the oscillatory difference between the two phase factors $\exp[{\mathrm{i}}k_{1}r_{j}]$ and $\exp[{\mathrm{i}}k_{2}r_{j}]$, leading to redundancy in the approximation space. Of course our method is designed with the high frequency regime in mind, and this redundancy is a low-frequency issue. But using slightly different grading parameters for the two meshes was found to be a simple and effective way to maintain stability even for relatively low frequencies.

4.2 Numerical experiments

In §4.3–§4.8 below, we present numerical results demonstrating the accuracy and efficiency of our HNA BEM for a range of scattering scenarios. As we shall see, our HNA BEM achieves an order of magnitude improvement in accuracy over the purely asymptotic PGOH approach, using only a modest (and frequency-independent) number of degrees of freedom.

To quantify the accuracy of the PGOH and HNA BEM solutions we compare them to reference solutions obtained by solving (11) using a conventional $hp$-BEM with a large number of degrees of freedom. Numerical experiments suggest that this reference solution achieves at least $0.01\%$ accuracy for all the problems considered here; for more details see [29]. We shall present results both for the boundary solution $\mathbf{v}=(u_{1},\partial{u_{1}}/\partial{\mathbf{n}})^{T}$ and the far-field pattern $F$ (see (10)).

Our algorithm can be applied to any convex polygonal scatterer, and to illustrate this we consider three different examples: an equilateral triangle, a square, and a regular hexagon. In each case, we take the sides of the polygon to be of length $2\pi$ so that the number of exterior (interior) wavelengths $\lambda_{1}$ ($\lambda_{2}$) around $\Gamma$ is equal to $n_{s}k_{1}$ ($n_{s}k_{2}=n_{s}k_{1}\mu$). We shall report results for $k_{1}$ ranging from $5$ to $160$, the highest value corresponding to $480$, $640$ and $960$ exterior wavelengths around $\Gamma$ for the triangle, square, and hexagon, respectively. (The reason we do not present results for even higher wavenumbers is that computing the reference solution becomes prohibitively expensive.) For comparison with conventional methods (such as the standard $hp$-BEM used to compute our reference solution) we introduce the notation

[TABLE]

to denote the average number of degrees of freedom (DOF) per wavelength around $\Gamma$ used in the approximation of each of the components of the boundary solution $\mathbf{v}$, relative to the exterior and interior wavelengths respectively. A rule of thumb for scattering problems (see, e.g., [43, 44]) is that to achieve acceptable “engineering accuracy” using a conventional method, one generally requires 6 to 10 degrees of freedom per wavelength. Precisely what is meant by “engineering accuracy” depends on the application at hand, but in [43, 44] the target accuracy appears to be roughly $10\%$. As we shall see, at high frequencies our HNA BEM can achieve better than $1\%$ accuracy with far less than one degree of freedom per wavelength.

For the equilateral triangle we shall consider scattering for the five different incident directions $\mathbf{d}_{1}$, $\mathbf{d}_{2}$, $\mathbf{d}_{3}$, $\mathbf{d}_{4}$ and $\mathbf{d}_{5}$ shown in Figure 7(a), which correspond to taking

[TABLE]

respectively. We also study the performance of our method under changes in the refractive index $\mu$, and the constant $\alpha$ in (6). Our work is motivated in part by applications in atmospheric physics, specifically the scattering of electromagnetic radiation by ice crystals in cirrus clouds [3]. In this context the relevant choices of $\alpha$ are $\alpha=1$ and $\alpha=(k_{1}/k_{2})^{2}$, as already discussed after (7). The refractive index $\mu$ of ice has been determined through numerous experiments (see, e.g., [54]) and is known to be complex-valued and highly frequency-dependent, with $\mu$ (the contrast) ranging between $1$ and $2$, and ${\rm Im}\left[\mu\right]$ (the absorption) between [math] and $1$. For most of our examples we take $\mu=1.5$. But we also present results for $\mu=0.66$ (the approximate reciprocal of $1.5$), to demonstrate that our method is also effective for the case $0<\mu<1$. We consider a range of ${\rm Im}\left[\mu\right]$ from [math] up to $0.0125$; results for larger absorptions (when our method achieves accuracies even better than those reported here) can be found in [29].

4.3 High frequency performance

In this section we investigate how the accuracy of our method depends on the wavenumber $k_{1}$. We consider the configuration from Figure 1, namely scattering by an equilateral triangle with $\mathbf{d}^{i}=\mathbf{d}_{1}$, $\alpha=1$ and $\mu=1.5+0.003125{\mathrm{i}}$, for wavenumbers $k_{1}=5,10,20,40,80,160$. The parameter choices (31)–(35) result in a total number of degrees of freedom $N=416$ in this case, independent of $k_{1}$.

Figure 8 shows the resulting error in our HNA BEM approximation to the boundary solution $\mathbf{v}=(u_{1},\partial{u_{1}}/\partial{\mathbf{n}})^{T}$ and the far-field pattern $F$. For a comparison with a purely asymptotic method we also show the corresponding errors for the GO approximation, and the resulting PGOH approximation to the far-field pattern (obtained by substituting the GO approximation for $\mathbf{v}$ into (10)). For a comparison with a purely numerical method we also show the corresponding errors for a “conventional” BEM, namely the same $hp$-BEM used to generate our reference solution, but with far fewer degrees of freedom; specifically, for the $hp$-BEM we use approximately the same ($k_{1}$-independent) number of degrees of freedom (in this case $N=456$) as we use in the HNA BEM ($N=416$). The corresponding values of $\mathrm{\#DOFper}\lambda_{2}$ for the two methods, along with the 2-norm condition numbers COND of the associated linear systems, are presented in Table 1.

For all values of $k_{1}$ considered here, the HNA method is accurate to within at least $5\%$ both for the computation of the boundary solution and the far-field pattern, and provides an order of magnitude improvement over GO/PGOH. Moreover, as is true for GO/PGOH, the accuracy of the HNA method increases as $k_{1}$ increases, despite the fact that $N$ remains fixed. For the largest wavenumber $k_{1}=160$, the HNA method achieves an accuracy of approximately $0.1\%$ using fewer than $0.3$ degrees of freedom per wavelength. In contrast, the error for the conventional $hp$-BEM (again with fixed $N$), despite starting below that of the HNA method for $k_{1}=5$, rapidly increases with increasing $k_{1}$, with the accuracy cross-over between $hp$-BEM and HNA BEM occurring between $k_{1}=10$ and $k_{1}=20$. The $hp$-BEM loses accuracy completely (relative error $>100\%$) for $k_{1}\geq 40$. (Of course, accuracy could be regained for the $hp$-BEM for $k_{1}\geq 40$ by increasing the number of DOF, but the point of Figure 8 is to compare the performance of the conventional $hp$-BEM and our HNA BEM when they are allocated a fixed (and roughly equal) number of DOFs.)

Table 1 reveals that, for the examples considered, the condition number for the HNA BEM broadly decreases as frequency increases, in contrast to the conventional $hp$-BEM. This is consistent with the fact that as $k_{1}$ increases our HNA approximation space becomes increasingly well-suited to approximating the boundary solution. The conditioning presents no issues for the direct linear solver used in the HNA BEM, as the error plots in Figure 8 demonstrate.

4.4 Varying the incident angle

We now investigate how the accuracy of the method depends on the incident wave direction. As in the last section we fix $\alpha=1$ and $\mu=1.5+0.003125{\mathrm{i}}$, but now consider the five evenly spaced incident angles $\mathbf{d}_{1},\ldots,\mathbf{d}_{5}$ defined by (36) (see Figure 7(a)). We already presented field plots for the case $\mathbf{d}^{i}=\mathbf{d}_{1}$ in Figure 1. Figure 7(b) shows an analogous field plot for the case $\mathbf{d}^{i}=\mathbf{d}_{3}$, which corresponds to grazing incidence along $\Gamma_{1}$. The number of degrees of freedom $N$ used in the approximation space for each of these examples ranges from 400 and 424; exact values are detailed in Table 2. The variation between examples is due to the variation of the strength of beam boundary discontinuities, and whether or not extra mesh points are introduced to capture them (recall (29), (24) and the discussion at the end of §3.1).

In Figure 9 we report $L^{2}$ errors in the approximations to $u_{1}$ on $\Gamma$ for all five incident directions and a range of wavenumbers $k_{1}$. (Errors in $\partial{u_{1}}/\partial{\mathbf{n}}$ and $F$ follow very similar trends but are not reported here.) In all cases considered, the HNA BEM provides a significant improvement over GO. It is interesting to note, however, that the errors in both the GO and HNA approximations are not uniform across incident angle; rather one observes a peak in the error curves at the grazing incidence case $\mathbf{d}^{i}=\mathbf{d}_{3}$ (i.e. $\theta_{i}=\pi/3$). We hypothesize that the larger error in the GO at grazing incidence is due to the fact that the high frequency asymptotic behaviour in this case involves a particularly prominent diffracted wave along $\Gamma_{1}$, which is not captured by the GO approximation. (By contrast, the HNA approximation does capture this diffracted wave.)

We hypothesize that the larger error in the HNA approximation at grazing incidence is due to a particularly prominent lateral wave contribution on $\Gamma_{3}$, which isn’t captured by our HNA approximation space. As was suggested in the caption to Figure 3, it is plausible that grazing incidence might lead to particularly large lateral wave contributions, because of the following argument: while the lateral wave usually compensates for the phase mismatch between the exterior and interior diffracted waves, at grazing incidence one of the lateral waves (in this case the one associated with $\Gamma_{1}$, since for $\mathbf{d}^{i}=\mathbf{d}_{3}$ the direction of the incident wave is parallel to that side) has to match up to the incident wave; since the incident wave has lower asymptotic order than the diffracted wave, one might therefore expect the lateral wave to be more prominent in this case.

Lateral waves only propagate inside $\Omega_{2}$ when $\mu>1$, which is the case for the examples considered up to now. When $0<\mu<1$ the lateral waves propagate away from the scatterer in $\Omega_{1}$, and therefore should not affect the boundary solution. Hence we can test our hypothesis about the role that lateral waves play in the deterioration of accuracy at grazing incidence for the case $\mu>1$ by studying the case $0<\mu<1$ and observing whether the same deterioration in accuracy occurs. This we do in the next section.

4.5 Varying the contrast ($\mu$)

In this section we study the performance of the HNA method in the case $0<\mu<1$. Specifically, fixing $\alpha=1$, we compare the results presented in Figure 9 for the case $\mu=1.5+0.003125{\mathrm{i}}$, with those presented in Figure 10 for $\mu=0.66+0.003125{\mathrm{i}}$ (the real part of the latter being approximately the reciprocal of the real part of the former). Based on these plots we make a number of observations. First, the accuracy of the GO approximation is roughly the same for the two refractive indices. Second, as was the case for $\mu=1.5$, the GO approximation for $\mu=0.66$ exhibits its largest errors at grazing incidence - again we ascribe this to a particularly prominent diffracted field along the grazing side in this case. Third, the HNA approximation consistently produces smaller errors for $\mu=0.66$ than for $\mu=1.5$. Fourth, the HNA error for $\mu=0.66$ is not obviously worse at grazing incidence than for other incident directions, in contrast to the $\mu=1.5$ case. To make a direct comparison easier, in Figure 11(a) we plot the errors for the two cases ($\mu=1.5$ and $\mu=0.66$) at the grazing incidence angle $\theta_{i}=\pi/3$ ($\mathbf{d}^{i}=\mathbf{d}_{3}$), for a range of values of $k_{1}$. The GO errors are virtually identical for the two cases, but the HNA errors are significantly smaller for $\mu=0.66$ than for $\mu=1.5$. These observations all support our hypothesis that the deterioration in accuracy of our HNA method near grazing incidence for $\mu>1$ is due to the presence of prominent lateral wave contributions, which are not present for $\mu<1$.

To test this hypothesis further, we consider the analogous experiment for a square scatterer. Aligning the square with the Cartesian axes, the grazing incidence case is now $\mathbf{d}^{i}=\mathbf{d}_{1}$ ($\theta_{i}=\pi/2$). Field plots for $u_{d}$ in this configuration for the two cases ($\mu=1.5$ and $\mu=0.66$) can be found in Figure 12. The corresponding boundary errors are presented in Figure 11(b). As for the triangle, the HNA error is significantly smaller for $\mu=0.66$ than for $\mu=1.5$, and again we hypothesize that this is due to the fact that for $\mu=1.5$ (and not for $\mu=0.66$) prominent lateral waves are generated by the two vertical sides, which impinge on the bottom side and limit the accuracy of our HNA approximation. We note also that this example allows us to rule out the possibility that increased error at grazing incidence might be due to the GO1/GO2 choice discussed in §3.1, because for this square configuration the propagation and decay directions ($\mathbf{d}$ and $\mathbf{e}$ in the notation of §3.1) always coincide and point vertically up or down, so that GO1 and GO2 coincide for all GO beams.

4.6 Varying $\alpha$

In principle, our HNA method can be evaluated for any value of the boundary condition jump parameter $\alpha\in\mathbb{C}\setminus\{0\}$. In practice, as was discussed after (7), the choices of $\alpha$ relevant to the dielectric electromagnetic scattering problem are $\alpha=1$ and $\alpha=(k_{1}/k_{2})^{2}$. We have already presented results for the former case; in this section we present results for the latter case. Specifically, we consider scattering by the triangle with $\mu=1.5+0.00625{\mathrm{i}}$ and $\mathbf{d}^{i}=\mathbf{d}_{1}$ ($\theta_{i}=\pi/2$). Figure 13 shows that for this configuration the performance of the method in the approximation of both $u_{1}$ and $\partial u_{1}/\partial\mathbf{n}$ is very similar for the two cases $\alpha=1$ and $\alpha=(k_{1}/k_{2})^{2}$.

4.7 Varying the absorption (${\rm Im}\left[\mu\right]$)

In this section we study the effect of varying the absorption ${\rm Im}\left[\mu\right]$. Specifically, we consider scattering by the triangle with $\mathbf{d}^{i}=\mathbf{d}_{5}$ ($\theta_{i}=\pi/6$), $\alpha=1$, $\mu=1.5$ and

[TABLE]

In Figure 14 we present the resulting error plots for the GO and HNA approximations. For all values of ${\rm Im}\left[\mu\right]$ considered (even ${\rm Im}\left[\mu\right]=0$) the HNA method provides an improvement in accuracy over the GO method and (once $k_{1}$ is sufficiently large) the HNA error decreases as $k_{1}$ increases. For fixed $k_{1}$ the HNA error decreases as the absorption ${\rm Im}\left[\mu\right]$ increases; for the worst case (${\rm Im}\left[\mu\right]=0$) using the HNA approximation instead of GO alone reduces the error at $k_{1}=10$ from $21.2\%$ to $2.8\%$ and at $k_{1}=160$ from $9.1\%$ to $4.0\%$, whereas for the best case (${\rm Im}\left[\mu\right]=0.0125$) the error at $k_{1}=10$ is reduced from $19.8\%$ to $2.0\%$ and at $k_{1}=160$ from $6.7\%$ to $0.17\%$. This dependence of accuracy on absorption can be explained in terms of the asymptotic theory. Our HNA approximation space neglects the contribution of the higher-order asymptotic terms such as lateral waves and reflected-diffrated waves (see Figure 3), and for ${\rm Im}\left[\mu\right]>0$ these higher-order terms are attenuated as they propagate across $\Omega_{2}$ at a rate proportional to ${\rm Im}\left[\mu\right]$. Hence the larger ${\rm Im}\left[\mu\right]$ is, the smaller the contribution the neglected terms make to the boundary solution $\mathbf{v}$, and, as a result, the more accurate our HNA BEM is.

4.8 Scattering by a hexagon

We have already presented results for scattering by an equilateral triangle and a square. In this section we consider scattering by a regular hexagon. This can be viewed as a simple two-dimensional model of the motivating application we mentioned in §1, namely the scattering of light by atmospheric ice crystals, which typically take the form of hexagonal columns or their aggregates (see [3]). Here we choose a refractive index $\mu=1.39+0.00667{\mathrm{i}}$ (which corresponds to the refractive index of ice for free-space wavelength 3.73$\mu m$ [53]), jump parameter $\alpha=1$ (corresponding to the out-of-plane electric field polarisation) and incident direction $\mathbf{d}^{i}=(1/\sqrt{13})(3,-2)$ (incident angle $\theta_{i}=\tan^{-1}(2/3)$). The total field and far-field pattern for this configuration are plotted in Figure 15.

The corresponding errors for $k_{1}=10,20,40,80,160$ are reported in Figure 16. The number of degrees of freedom $N$ in the HNA method decreases from 1008 at $k_{1}=10$ to 912 at $k_{1}=160$. (As was remarked in §4.4, our method typically uses fewer degrees of freedom at larger wavenumbers because beam boundary discontinuities decrease in strength as the wavenumber increases.) For all wavenumbers considered, the HNA method provides a significant improvement in accuracy over GO (PGOH in the far field), and increases in accuracy as $k_{1}$ increases, despite the fact that the number of degrees of freedom $N$ remains fixed or decreases. Specifically, for the approximation of the far-field pattern at the highest wavenumber $k_{1}=160$, which corresponds to $960$ exterior wavelengths (approximately $1334$ interior wavelengths) around the scatterer boundary, the PGOH approximation achieves $1.9\%$ error, while the HNA BEM achieves $0.049\%$ error with just $N=912$, which corresponds to $\mathrm{\#DOFper}\lambda_{2}\approx 0.34$.

5 Discussion

In this paper we have described an HNA BEM for high frequency scattering by penetrable (dielectric) convex polygons. We have demonstrated, by a range of numerical experiments, that our method can provide an order of magnitude improvement in accuracy over the purely asymptotic GO (PGOH in the far field) approach. Moreover, our method can achieve a fixed accuracy of approximation using a relatively small, frequency-independent number of degrees of freedom $N$. By contrast, conventional BEM approaches require $N$ to grow at least linearly with increasing frequency in order to maintain accuracy. The extremely efficient deployment of degrees of freedom in our HNA BEM is due to the fact that it uses an approximation space built from oscillatory basis functions carefully chosen to capture certain components of the high frequency solution behaviour.

We conclude the paper with a discussion of some possible (but non-trivial) modifications that might enhance the performance of our method still further. First we recall that (as was explained in detail in §3) our HNA method, while more accurate than GO alone, is not fully error-controllable (in the sense of the HNA methods analysed rigorously in [17, 35, 16, 33, 34]), because its accuracy is limited by our neglect of higher-order asymptotic effects such as lateral waves and diffracted-reflected waves (see Figure 3). Full HNA error-controllability seems infeasible for this problem because there are infinitely many different phases to consider, despite the fact that the scatterer is convex. But, in principle, the accuracy of our method could be improved, at the expense of additional degrees of freedom, by incorporating new oscillatory basis functions capable of capturing some of these higher-order effects. Concretely, for the lateral wave contributions one could include plane-wave basis functions of the form $f(\mathbf{x}){\mathrm{e}}^{{\mathrm{i}}k_{2}\mathbf{d}^{l}\cdot\mathbf{x}}$, where $f$ is a piecewise polynomial and $\mathbf{d}^{l}$ is the propagation direction of the lateral wave in question (these directions can be determined solely from the refractive index and the corner angles of the polygon - for further discussion of these issues see [31, §3.2.1]). For the diffracted-reflected waves one could include basis functions of the form $f(\mathbf{x}){\mathrm{e}}^{{\mathrm{i}}k_{2}r_{j}^{\prime}(\mathbf{x})}$, where $f$ is a piecewise polynomial and $r_{j}^{\prime}(\mathbf{x})$ is the distance between $\mathbf{x}$ and an appropriate reflected image $\mathbf{P}_{j}^{\prime}$ of the corner $\mathbf{P}_{j}$ from which the original diffracted ray field emanated. (Such an image point $\mathbf{P}_{j}^{\prime}$ is shown in Figure 3(b).) We note that similar basis functions are required to capture multiple scattering effects in HNA BEMs for nonconvex impenetrable polygons - see, e.g., the discussion in [16, §8].

Second, we suspect that our rather simplistic treatment of GO beam boundaries may not be optimal. It is well-known that, for scattering by impenetrable wedges, the smooth (but rapidly-varying) solution behaviour near the “shadow boundaries” across which GO components (incident or reflected) switch on/off is governed by the Fresnel integral and related functions (see, e.g., [48, 8]). In the context of HNA methods for sound-soft nonconvex polygons this (frequency-dependent) shadow boundary behaviour can be captured either by including the appropriate special functions in the HNA approximation space (as in [16]) or by cutting off the GO component sharply across the shadow boundary and refining the mesh for the associated diffracted component towards it, to capture the rapid variation that compensates for the GO discontinuity (as in [33]). Our approach in the current paper, in the context of penetrable convex polygons, is a rather crude implementation of the latter approach. We expect that certain minor (yet technical) modifications to our current practice might lead to slightly improved accuracy:

• •

Following [31], we take the beam boundaries of a GO beam $a{\mathrm{e}}^{{\mathrm{i}}(D\mathbf{d}+{\mathrm{i}}E\mathbf{e})\cdot\mathbf{x}}$ to be parallel to the associated propagation vector $\mathbf{d}$. When $E>0$ (which will be the case generically if ${\rm Im}\left[\mu\right]>0$) this choice may not be optimal - in simpler “knife-edge” diffraction problems involving inhomogeneous incident waves the “correct” location of shadow boundaries (defined in an appropriate sense) is known to be shifted in a rather non-intuitive way (see [4] and the discussion in [31, Remark 3.1]). A possible modification to our current algorithm would be to allow an absorption-dependent shift in the beam boundary locations.

• •

For each beam boundary discontinuity classified as “strong” (in the sense of (24)), we insert a single new point in the mesh for the associated diffracted wave (if one exists). The rigorous analysis in [33] for sound-soft nonconvex polygons suggests that the rapid solution variation near the beam boundary might be captured more accurately by introducing geometric mesh grading towards the beam boundary.

To our knowledge, the canonical penetrable wedge scattering problem has not been analysed in sufficient detail to allow us to make theoretically-justifiable judgements about the optimal strategy in relation to the points raised above. Hence we leave such considerations to future work.

Funding/acknowledgements

This work was supported by EPSRC (EP/K000012/1 to SL, EP/P511262/1 to DPH, PhD studentship to SPG) and the Met Office (PhD CASE award to SPG). We would like to thank Anthony Baran for helpful comments and suggestions.