# Integral equation methods for acoustic scattering by fractals

António Caetano, Simon Chandler-Wilde, Xavier Claeys, Andrew Gibbs, David Hewett, Andrea Moiola

► To cite this version:

António Caetano, Simon Chandler-Wilde, Xavier Claeys, Andrew Gibbs, David Hewett, et al. Integral equation methods for acoustic scattering by fractals. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2025, 481 (2306), 10.1098/rspa.2023.0650. hal-05497675

HAL Id: hal-05497675

https://hal.science/hal-05497675v1

Submitted on 6 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Integral equation methods for acoustic scattering by fractals

A. M. Caetano*, S. N. Chandler-Wilde†, X. Claeys‡,

A. Gibbs§, D. P. Hewett§ and A. Moiola¶

## Abstract

We study sound-soft time-harmonic acoustic scattering by general scatterers, including fractal scatterers, in 2D and 3D space. For an arbitrary compact scatterer $\Gamma$ we reformulate the Dirichlet boundary value problem for the Helmholtz equation as a first kind integral equation (IE) on $\Gamma$ involving the Newton potential. The IE is well-posed, except possibly at a countable set of frequencies, and reduces to existing single-layer boundary IEs when $\Gamma$ is the boundary of a bounded Lipschitz open set, a screen, or a multi-screen. When $\Gamma$ is uniformly of $d$-dimensional Hausdorff dimension in a sense we make precise (a $d$-set), the operator in our equation is an integral operator on $\Gamma$ with respect to $d$-dimensional Hausdorff measure, with kernel the Helmholtz fundamental solution, and we propose a piecewise-constant Galerkin discretization of the IE, which converges in the limit of vanishing mesh width. When $\Gamma$ is the fractal attractor of an iterated function system of contracting similarities we prove convergence rates under assumptions on $\Gamma$ and the IE solution, and describe a fully discrete implementation using recently proposed quadrature rules for singular integrals on fractals. We present numerical results for a range of examples and make our software available as a Julia code.

## 1 Introduction

This paper, prepared in large part during a recent Isaac Newton Institute programme on multiple wave scattering, is concerned with the classical prob-

<footnote>*Center for R&D in Mathematics and Applications, Departamento de Matemática, Universidade de Aveiro, Aveiro, Portugal</footnote>

<footnote>†Department of Mathematics and Statistics, University of Reading, Reading, United Kingdom</footnote>

<footnote>‡Laboratoire Jacques-Louis Lions, Sorbonne Université, Paris, France</footnote>

<footnote>§Department of Mathematics, University College London, London, United Kingdom</footnote>

<footnote>¶Dipartimento di Matematica “F. Casorati”, Università degli studi di Pavia, Pavia, Italy</footnote>
lem of scattering of time-harmonic acoustic waves in $\mathbb{R}^n$, $n = 2$ or 3, by a scatterer $\Gamma$ (assumed to be a compact subset of $\mathbb{R}^n$) that may have multiple components or other complicated geometrical features. We consider the sound-soft case, where the total field $u^t$ vanishes on $\Gamma$, and satisfies the Helmholtz equation $(\Delta + k^2)u^t = 0$ for some wavenumber $k > 0$ in the open set $\Omega := \mathbb{R}^n \setminus \Gamma$. Our focus is on integral equation (IE) formulations of this scattering problem and their numerical solution<fnref n="1" />. Our particular interest is in scattering by fractals, which provide a model for the multiscale roughness of many natural and man-made scatterers.

In a sequence of recent papers [9, 11, 5] we studied scattering by fractal planar screens (i.e. scattering in $\mathbb{R}^n$ by fractal subsets of $\mathbb{R}^{n-1}$). In the current paper we show how the results of [9, 11, 5] can be generalised to scattering by fractal subsets of $\mathbb{R}^n$ that are not contained in a hyperplane. This is a significant novelty compared to previous contributions, greatly extending the class of scatterers to which our results apply. This generalisation complicates the analysis, because when the scatterer is a planar screen as in [9, 11, 5], irrespective of its smoothness, the scattering problem can be written as a coercive (sign-definite) variational problem by the results in [8, 19, 20], while in the general case considered here the integral operators involved are compact perturbations of coercive ones and we have to resort to Fredholm theory. This generalisation also forces us to use a novel Sobolev space setting for the integral equation. On the other hand, for the case where the scatterer is a planar screen, the analysis carried out here is in certain respects a simplification of that in [9, 11, 5], as the trace operator from $\mathbb{R}^n$ to a hyperplane (see Lemma 3.11) does not play a role.

The contributions of the paper are several. Firstly, we write down in Theorem 3.4 a novel IE formulation of the problem that applies for any compact $\Gamma \subset \mathbb{R}^n$. The scattered field is sought as an acoustic Newton potential $u = A\phi$, with an unknown density $\phi$ on $\Gamma$ satisfying a first-kind IE

$$A\phi = g, \qquad (1)$$

for some data $g$ depending on the incident wave. The operator $A$, defined in (16) below, is a generalisation of the classical single-layer boundary integral operator, in cases where this is well-defined. Importantly, we prove in Lemma 3.3 that $A$ is a compact perturbation of a coercive operator, so that all Galerkin solution methods are convergent, provided $A$ is also in- jective. Using this result, in Theorem 3.4 we provide an IE-based proof of well-posedness for scattering by a general compact $\Gamma$, generalising existing IE-based proofs for cases where $\Omega$ is Lipschitz or smoother (e.g., [29, Thm 9.11]).

<footnote n="1">¹We note that our methods and results apply, with obvious modifications, to the analogous (yet simpler) problem in potential theory, in which the Helmholtz equation is replaced by the Laplace equation.</footnote>
We focus mostly on the particular case where $\Gamma$ is a $d$-set (definition (2) below), which means (roughly speaking) that $\Gamma$ is uniformly of Hausdorff dimension $d$, for some integer or fractional $d \in (0, n]$. If $\Gamma$ is a $d$-set and $d \le n-2$, then, as we explain in Remark 3.5, the scatterer is invisible to incident waves. To focus on cases where $u \ne 0$ we restrict our study to the range $n-2 < d \le n$.

A range of examples with $n-2 < d \le n$, relevant to our later discussions and computations, is pictured in Figures 1 and 2. Figure 1 shows examples in 2D space ($n=2$), namely: (a) $\Gamma = \bar{D}$ is the closure of a bounded Lipschitz open set $D$ ($d=n=2$); (b) $\Gamma = \partial D$ is the boundary of the same set ($d=n-1=1$); (c) $\Gamma = [0,1] \times \{0\}$ is a line segment ($d=n-1=1$); (d) $\Gamma = [-1,1] \times \{0\} \cup \{0\} \times [-1,1]$ is the cross formed by two line segments, an example of a multi-screen in the sense of [13] (all such multi-screens are $d$-sets with $d=n-1$); (e) $\Gamma = C \times \{0\}$, where $C \subset [0,1]$ is the classical middle-third Cantor set ($d = \log(2)/\log(3) \approx 0.6309$); (f) $\Gamma$ is the Koch curve ($d = \log(4)/\log(3) \approx 1.262$); (g) $\Gamma$ is the closure of the Koch snowflake domain ($d=2$).

Examples (c) and (e)-(g) in Figure 1 are all fixed points (attractors) of an iterated function system (IFS) satisfying the standard open set condition (OSC) (we recall these definitions in §2.1). As we recall in §2.1, every such IFS attractor is a $d$-set, with $d$ its fractal (Hausdorff) dimension. Figure 2 shows examples in 3D space ($n=3$) where $\Gamma$ is an IFS attractor that is a Sierpinski tetrahedron, with $d=2$ and $d = \log 4 / \log(8/3) \approx 1.41$. (We show numerical simulations for scattering by these shapes in §5.) These examples make clear that our results include multiple scattering cases where $\Gamma$ has a complicated geometry and/or multiple components. Indeed, the IFS attractor examples in Figure 1(e) and Figure 2(b) are both fractal cases where the IFS is disjoint (as defined in §2.1) so that $\Gamma$ is totally disconnected and has uncountably many components!

A key result, proved in Theorem 3.16, important both theoretically and computationally, is that in the $d$-set case we can interpret the Newton potential $A\phi$ as an integral with respect to $H^d$, the $d$-dimensional Hausdorff measure. Similarly, we show in Theorem 3.16 that the operator $A$ in (1) can also be written equivalently as an integral operator $A$ with respect to the $H^d$ measure.

In certain special cases our formulation coincides with previously studied IE formulations. If $d=n$ (e.g., Figure 1(a), (g)), so that $H^d$ is $n$-dimensional Lebesgue measure, $A$ is a volume integral operator and (1) is equivalent to a volume IE on $\Gamma$. (Note however that, where $\Omega_+$ is the unbounded component of $\Omega$, the solution of (1) is always supported in $\partial\Omega_+ \subset \partial\Gamma$; see Remark 3.7.) If $d=n-1$ and $\Gamma$ is the whole or part of the boundary of a bounded Lipschitz open set (e.g., Figure 1(b), (c)), or is a multi-screen in the sense of [13] (e.g., Figure 1(d)), then $H^d$ is standard surface measure (e.g., [15, Theorem 3.8]), $A$ is a surface integral operator, specifically an acoustic
<figure><img src="image_2.png" /><figcaption>Figure 1: Examples of *d*-sets in two-dimensional space (*n* = 2), with: a) *d* = 2; b) *d* = 1; c) *d* = 1; d) *d* = 1; e) *d* = log(2)/log(3) ≈ 0.63; f) *d* = log(4)/log(3) ≈ 1.26; g) *d* = 2. For details see text of §1.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Sierpinski tetrahedron $d$-sets in 3D space, attractors of the IFS (5) for $\rho = 1/2$ ($d = 2$) and $\rho = 3/8$ ($d = \log 4 / \log(8/3) \approx 1.41$).</figcaption></figure>

single-layer boundary integral operator, and (1) is equivalent to a standard first kind boundary IE (see Remark 3.17). Finally, in the case when $\Gamma$ is a subset of a hyperplane (e.g., Figure 1(c), (e)), our formulation reduces to cases studied recently in [10, 9, 11, 5, 4]. In particular, as already noted above, our results and methods build on those in [5, 4], where we study scattering by fractal planar screens that are $d$-sets with $n - 2 < d \le n - 1$.

In the $d$-set case we show that the integral operator $\mathcal{A}$ is a continuous mapping on a scale of Sobolev spaces on $\Gamma$ (Proposition 3.18), a first step in a regularity theory for solutions of (1). In the $d$-set case we are also able to propose (in §4, Eqn. (43) in particular) a piecewise-constant Galerkin IE method (IEM) for the numerical solution of (1), which we prove (in Theorem 4.3) is convergent as the mesh size $h$ tends to zero. Moreover, the entries of the matrix and right-hand-side of the linear system defining the Galerkin solution are given explicitly as double and single integrals, respectively, with respect to $H^d$ measure.

Our IEM is familiar in the case where $\Gamma$ is the whole or part of the boundary of a bounded Lipschitz open set, or is a multi-screen in the sense of [13]; our IEM is then a standard Galerkin boundary element method (BEM) with a piecewise-constant approximation space. In the case where $\Gamma$ is a $d$-set that is a planar screen our IEM coincides with that of [5]; indeed, our linear system is identical to that in [5, Eqn (55)].

Our strongest results (see §4.1) are for the special $d$-set case where $\Gamma$ is the attractor of an IFS satisfying the OSC (e.g., Figure 1(c), (e)-(g), Figure 2). In this case, under appropriate assumptions on $\Gamma$, we are able to prove (in Theorems 4.4 and 4.5) convergence rates for our IEM. The rate of
convergence depends on the regularity of the solution $\phi$. In Hypothesis 3.21 and Remark 4.6 we introduce a hypothesis relating to this regularity and detail the resulting convergence rates. In this case we also propose in §4.2 a fully discrete implementation, evaluating our Hausdorff single and double integrals using recently proposed quadrature methods for singular integrals on IFS attractors satisfying the OSC [18, 17].

In §5 we show computations using our fully discrete Galerkin IEM, solving (1) and computing the scattered field $u = A\phi$ for scatterers $\Gamma$ including the examples in Figure 2. This section includes numerical experiments exploring the convergence of our method. These suggest that our regularity hypothesis (Hypothesis 3.21) is true for many of the examples we study, and that the conditions of Theorem 4.4, guaranteeing the validity of our convergence rate analysis, may be satisfied generally whenever $\Gamma$ is an IFS attractor satisfying the OSC (establishing this, and proving some version of our regularity hypothesis, are both open problems).

We end this introduction by noting that an alternative approach to the simulation of scattering by fractals is to approximate the fractal by a smoother “prefractal” scatterer and apply a more conventional numerical method to the resulting approximate scattering problem. This was the approach taken in [11, 2], and in the earlier work in [22] and [31] for Laplace and elasticity problems. An achievement of [11, 2] was to prove convergence (without rates) of conventional BEMs for prefractal approximations of fractal planar screen problems, via Mosco convergence techniques. In principle, a similar analysis could be carried out for the problems under consideration in the current paper. However, we do not pursue this here.

## 2 Preliminaries

In this section we set some notation/terminology, and briefly review a number of known results that we will use later in the paper. Further details can be found in the references provided.

Throughout, for $n \in \mathbb{N}$ and $E \subset \mathbb{R}^n$, $\bar{E}$, $\partial E$, and $E^\circ := \bar{E} \setminus \partial E$ denote the closure, boundary, and interior of $E$ with respect to the standard Euclidean metric on $\mathbb{R}^n$, and $E^c := \mathbb{R}^n \setminus E$ its complement in $\mathbb{R}^n$. In the case that $E$ is measurable, $m(E)$ denotes its $n$-dimensional Lebesgue measure. $B_r(x) \subset \mathbb{R}^n$ denotes the closed ball of radius $r$ centred on $x$.

### 2.1 Hausdorff measure and dimension, $d$-sets, IFS attractors

For $0 \le d \le n$ let $\mathcal{H}^d$ denote the Hausdorff $d$-measure on $\mathbb{R}^n$ and let $\dim_H(S) \in [0, n]$ denote the Hausdorff dimension of $S \subset \mathbb{R}^n$ (see, e.g., [16]). For convenience we adopt the normalisation of [15, Def. 2.1], so that $\mathcal{H}^d$ coincides with Lebesgue measure for $d = n$. As in [24, §1.1] and [36,
§3], given $0 < d \le n$, we say a closed set $\Gamma \subset \mathbb{R}^n$ is a $d$-set if there exist $c_2 > c_1 > 0$ such that

$$c_1 r^d \le \mathcal{H}^d(\Gamma \cap B_r(x)) \le c_2 r^d, \quad x \in \Gamma, \ 0 < r \le 1. \quad (2)$$

We note that $d$-sets are also termed Ahlfors $d$-regular or Ahlfors-David $d$-regular sets, e.g., [27, p. 92]. We note also that if $\Gamma$ is a $d$-set then $\dim_H(\Gamma) = d$.

By an *iterated function system of contracting similarities* (we abbreviate this whole phrase as *IFS*)<fnref n="2" /> we mean a collection $\{s_1, s_2, \dots, s_M\}$, for some $M \ge 2$, where, for each $m = 1, \dots, M$, $s_m : \mathbb{R}^n \to \mathbb{R}^n$, with $|s_m(x) - s_m(y)| = \rho_m|x - y|$, $x, y \in \mathbb{R}^n$, for some $\rho_m \in (0, 1)$. The attractor of the IFS is the unique non-empty compact set $\Gamma$ satisfying

$$\Gamma = s(\Gamma), \quad \text{where} \quad s(E) := \bigcup_{m=1}^{M} s_m(E), \quad E \subset \mathbb{R}^n. \quad (3)$$

We shall restrict our attention to OSC-IFSs, i.e. IFSs that satisfy the standard open set condition (OSC) [16, (9.11)], which implies that the attractor $\Gamma$ is a $d$-set (see, e.g., [36, Thm. 4.7]), where $d \in (0, n]$ is the unique solution of $\sum_{m=1}^{M} (\rho_m)^d = 1$. For a homogeneous OSC-IFS, where $\rho_m = \rho \in (0, 1)$ for $m = 1, \dots, M$, we have $d = \log(M)/\log(1/\rho)$. If an OSC-IFS is not homogeneous we say it is non-homogeneous. Returning to the general, not necessarily homogeneous case, the OSC also implies (again, see [36, Thm. 4.7]) that $\Gamma$ is self-similar, meaning that the sets

$$\Gamma_m := s_m(\Gamma), \quad m = 1, \dots, M, \quad (4)$$

satisfy $\mathcal{H}^d(\Gamma_m \cap \Gamma_{m'}) = 0$, $m \ne m'$, so that $\Gamma$ is decomposed by (3) into $M$ similar copies of itself whose pairwise intersections have Hausdorff measure zero. For many of our results we make the additional assumption that the sets $\Gamma_1, \dots, \Gamma_M$ are disjoint. If this holds we say that the IFS attractor $\Gamma$ is disjoint, the OSC is automatically satisfied (e.g., [5, Lem. 2.5]), and $d < n$ (e.g. [5, Lemma 2.6]). If $\Gamma$ is not disjoint we say it is non-disjoint.

The following construction makes clear that if $C$ is an IFS attractor in dimension $n-1$ then $C \times \{0\}$ is an IFS attractor in dimension $n$.

**Remark 2.1 (Lifting attractors to higher dimensions).** Suppose that $M \ge 2$ and that $S = \{s_1, s_2, \dots, s_M\}$ is an IFS on $\mathbb{R}^n$. For $m = 1, \dots, M$, define $\tilde{s}_m : \mathbb{R}^{n+1} \to \mathbb{R}^{n+1}$ by $\tilde{s}_m((x, t)) = (s_m(x), \rho_m t)$, for $x \in \mathbb{R}^n, t \in \mathbb{R}$. Then, for $m = 1, \dots, M$, $\tilde{s}_m$ is a contracting similarity with the same contraction factor $\rho_m$ as $s_m$, so that $\tilde{S} = \{\tilde{s}_1, \tilde{s}_2, \dots, \tilde{s}_M\}$ is an IFS on $\mathbb{R}^{n+1}$. Further, $\tilde{S}$ satisfies the OSC/is disjoint if the same holds for $S$. If $\Gamma$ is the

<footnote n="2">²A useful introduction to IFSs is [16, Chap. 9]; the website [32] gives many examples of IFSs and their attractors.</footnote>
attractor of S then $\tilde{\Gamma} := \Gamma \times \{0\} = \{(x, 0) : x \in \Gamma\}$ is the attractor of $\tilde{S}$, and if the OSC holds for S then $\Gamma$ and $\tilde{\Gamma}$ are both d-sets with the same value of d.

**Example 2.2 (Cantor set examples of IFS attractors).** Let $S = \{s_1, s_2\}$, where $s_m : \mathbb{R} \to \mathbb{R}$, $m = 1, 2$, are defined, for some $\rho \in (0, 1/2]$, by

$$s_1(t) = \rho t, \quad s_2(t) = 1 + \rho(t-1), \quad t \in \mathbb{R}.$$

Then S is a homogeneous IFS with attractor C that is the “middle-(1-2ρ)” Cantor set, given by

$$C = \bigcap_{n=0}^{\infty} C_n, \quad \text{where} \quad C_0 := [0, 1], \quad C_n := s(C_{n-1}), \quad n \in \mathbb{N},$$

and s is the mapping given by (3) on the set of subsets of $\mathbb{R}$. In the case $\rho = 1/2$, $C = C_n = [0, 1]$ for each n. If $0 < \rho < 1/2$ then $C_n$ is a union of $2^n$ disjoint closed intervals and $C_n$ is obtained from $C_{n-1}$ by removing the middle $(1-2\rho)$ from each of the intervals comprising $C_{n-1}$. This IFS satisfies the OSC (see [16, §9.2]), so that (see above) it is a d-set with $d = \dim_H(C) = \log(2)/\log(1/\rho)$. The attractor C is disjoint if $\rho < 1/2$. By Remark 2.1, $C \times \{0\} \subset \mathbb{R}^2$, shown for $\rho = 1/2$ and $\rho = 1/3$ in Figure 1(c) and (e), respectively, is also the attractor of an IFS, is a d-set with the same value of d, and is disjoint if $\rho < 1/2$.

**Example 2.3 (Koch curve).** The Koch curve $\Gamma \subset \mathbb{R}^2$, shown in Figure 1(f), is the attractor of the homogeneous IFS $\{s_1, s_2, s_3, s_4\}$, where the mappings $s_m : \mathbb{R}^2 \to \mathbb{R}^2$, $m \in \{1, \dots, 4\}$, are given by

$$s_1(x) = x/3, \qquad s_2(x) = Rx/3 + \begin{pmatrix} \frac{1}{3} \\ 0 \end{pmatrix}, \\ s_3(x) = R^{-1}x/3 + \begin{pmatrix} \frac{1}{2} \\ \frac{1}{2\sqrt{3}} \end{pmatrix}, \qquad s_4(x) = x/3 + \begin{pmatrix} \frac{2}{3} \\ 0 \end{pmatrix},$$

for $x \in \mathbb{R}^2$, where R is the (orthogonal) rotation matrix for rotation counterclockwise by angle $\pi/3$. S satisfies the OSC (e.g., [16, Ex. 9.5]) and so $\Gamma$ is a d-set with $d = \log(4)/\log(3) \approx 1.26$. By Remark 2.1, $\Gamma \times \{0\} \subset \mathbb{R}^3$ is also the attractor of a homogeneous IFS, and is a d-set with the same value of d.

**Example 2.4 (Koch snowflake).** The Koch snowflake $\Gamma \subset \mathbb{R}^2$, shown in Figure 1(g), is the attractor of a non-homogeneous IFS of 7 contracting similarities satisfying the OSC. It is non-disjoint and has dimension $d = 2$. For details see [17, §5.4].
**Example 2.5 (Sierpinski tetrahedron).** A Sierpinski tetrahedron $\Gamma \subset \mathbb{R}^3$ can be defined, for every $0 < \rho \le 1/2$, as the attractor of the IFS comprising the four contracting similarities

$$s_i(x) := x_i + \rho(x - x_i), \quad i = 1, \dots, 4, \qquad (5)$$

where the $x_i$ are the vertices of a unit tetrahedron, explicitly

$$x_1 = (0, 0, 0)^T, \quad x_2 = (1, 0, 0)^T, \quad x_3 = (1/2, \sqrt{3}/2, 0)^T, \\ x_4 = (1/2, 1/(2\sqrt{2}), \sqrt{5}/(2\sqrt{2}))^T.$$

$\Gamma$ is shown in Figure 2 for $\rho = 3/8$ and $\rho = 1/2$. It satisfies the OSC, has dimension $d = \log 4 / \log(1/\rho)$, and is disjoint for $0 < \rho < 1/2$ but not for $\rho = 1/2$.

## 2.2 Function spaces

In this section we briefly review some function space definitions and results that will be used in our scattering problem and its IE formulations. Our notation follows that of [29] and [5]. Throughout, our function spaces are spaces of complex-valued functions/distributions.

For $s \in \mathbb{R}$ we let $H^s(\mathbb{R}^n)$ denote the usual Bessel potential Sobolev space. For a non-empty open set $\Omega \subset \mathbb{R}^n$, where $C_0^\infty(\Omega)$ is the set of those $C^\infty$ functions that are compactly supported in $\Omega$, we define $\tilde{H}^s(\Omega) := \overline{C_0^\infty(\Omega)}^{H^s(\mathbb{R}^n)}$, a closed subspace of $H^s(\mathbb{R}^n)$. For a non-empty closed set $E \subset \mathbb{R}^n$ we denote by $H_E^s$ the set of all elements of $H^s(\mathbb{R}^n)$ whose support is contained in $E$, also a closed subspace of $H^s(\mathbb{R}^n)$. We recall that, for $s \in \mathbb{R}$, $H^{-s}(\mathbb{R}^n)$ is dual to $H^s(\mathbb{R}^n)$, with the duality pairing $\langle \cdot, \cdot \rangle_{H^{-s}(\mathbb{R}^n) \times H^s(\mathbb{R}^n)}$ extending the $L_2(\mathbb{R}^n)$ inner product,<fnref n="3" /> and that, with respect to this same pairing, $H_E^{-s}$ is dual to $\tilde{H}^s(E^c)^\perp$, the orthogonal complement of $\tilde{H}^s(E^c)$ in $H^s(\mathbb{R}^n)$, for any non-empty closed $E \subsetneq \mathbb{R}^n$ [10, Cor. 3.4].

For non-empty open sets $\Omega \subset \mathbb{R}^n$ we also work with the classical Sobolev space $W^1(\Omega)$, normed by $\|u\|_{W^1(\Omega)}^2 = \|u\|_{L_2(\Omega)}^2 + \|\nabla u\|_{L_2(\Omega)}^2$, the closed subspace $W_0^1(\Omega) = \overline{C_0^\infty(\Omega)}^{W^1(\Omega)}$, and their “local” versions $W^{1,\text{loc}}(\Omega)$ and $W^{1,\text{loc}}(\Omega)$, defined as the sets of measurable functions $v$ on $\Omega$ such that $\sigma|_{\Omega}v$ is in $W^1(\Omega)$ or $W_0^1(\Omega)$, respectively, for every $\sigma \in C_0^\infty(\mathbb{R}^n)$. Similarly, we define local versions $H^{1,\text{loc}}(\mathbb{R}^n)$ and $\tilde{H}^{1,\text{loc}}(\Omega)$ of $H^1(\mathbb{R}^n)$ and $\tilde{H}^1(\Omega)$. We note that $H^1(\mathbb{R}^n) = W^1(\mathbb{R}^n)$ (and hence $H^{1,\text{loc}}(\mathbb{R}^n) = W^{1,\text{loc}}(\mathbb{R}^n)$), and that $\tilde{H}^1(\Omega) = W_0^1(\Omega)$ (and hence $\tilde{H}^{1,\text{loc}}(\Omega) = W_0^{1,\text{loc}}(\Omega)$) for arbitrary non-empty open $\Omega \subset \mathbb{R}^n$, with the latter identification involving the restriction operator, with extension by zero as its inverse.

Finally, for compact $\Gamma \subset \mathbb{R}^n$ let $C_{0,\Gamma}^\infty$ denote the set of functions in $C_0^\infty(\mathbb{R}^n)$ that equal one in a neighbourhood of $\Gamma$.

<footnote n="3">³Note that all our distributions and dual spaces are anti-linear rather than linear to suit our complex Hilbert space setting.</footnote>
### 3 Scattering problem and integral equation formulations

Let $\Gamma \subset \mathbb{R}^n$ ($n = 2,3$) be non-empty and compact and let $k > 0$. We consider the time-harmonic acoustic scattering of an incident wave $u^i$ by $\Gamma$, a sound-soft obstacle. We assume that the incident wave $u^i$ is an element of $W^{1,\text{loc}}(\mathbb{R}^n) = H^{1,\text{loc}}(\mathbb{R}^n)$ satisfying the Helmholtz equation

$$ \Delta u + k^2 u = 0 \qquad (6) $$

in a distributional sense in some neighbourhood of $\Gamma$ (so that $u^i$ is $C^\infty$ in that neighbourhood by elliptic regularity, see, e.g., [14, Thm 6.3.1.3]); for instance, $u^i$ might be the plane wave $u^i(x) = e^{ik\vartheta \cdot x}$ for some $\vartheta \in \mathbb{R}^n$ with $|\vartheta| = 1$). Where $\Omega := \Gamma^c = \mathbb{R}^n \setminus \Gamma$, we seek a scattered field $u \in W^{1,\text{loc}}(\Omega)$ satisfying (6) in a distributional sense in $\Omega$ (so that $u \in C^\infty(\Omega)$ by elliptic regularity), the Sommerfeld radiation condition

$$ \frac{\partial u(x)}{\partial r} - iku(x) = o(r^{-(n-1)/2}), \quad r := |x| \to \infty, \text{ uniformly in } \hat{x} := x/|x|, \qquad (7) $$

and the boundary condition $u = -u^i$ on $\partial\Omega = \partial\Gamma$, enforced by requiring that the total field

$$ u^t := u + u^i \in W_0^{1,\text{loc}}(\Omega). \qquad (8) $$

Note that $\partial\Omega = \Gamma$ if and only if $\Gamma$ has empty interior.

This problem, which we will refer to as our scattering problem, is uniquely solvable in the case that $\Omega$ is connected (see, e.g. [12, §3]). Figure 1, with the exception of (b), and Figure 2 are all examples of such cases. But, to understand the well-posedness of our IE formulation, we also want to allow cases (such as Figure 1(b)) where $\Omega$ is not connected, in which case $\Omega = \Omega_+ \cup \Omega_-$, where $\Omega_\pm$ are disjoint open sets, with $\Omega_+$ the unbounded component of $\Omega$ and $\Omega_-$ a bounded open set. In such cases the above problem decouples into a uniquely-solvable scattering problem for $u|_{\Omega_+} \in W^{1,\text{loc}}(\Omega_+)$ and the homogeneous Dirichlet problem that $u^t|_{\Omega_-} \in W_0^1(\Omega_-)$ satisfies (6) in $\Omega_-$. Thus, if $\Omega$ is not connected, our scattering problem is uniquely solvable (with $u^t = 0$ and $u = -u^i$ in $\Omega_-$) if and only if $k^2$ is not a Dirichlet eigenvalue of $-\Delta$ in $\Omega_-$. We will frequently assume that $k$ is not one of these exceptional values, making the following assumption.

**Assumption 3.1.** The only $v \in W_0^{1,\text{loc}}(\Omega)$ satisfying (6) in $\Omega := \Gamma^c$ and (7) is $v = 0$.

**Remark 3.2.** We emphasise that Assumption 3.1 holds for all $k > 0$ if $\Omega$ is connected, and that if $\Omega$ is not connected, in which case $\Omega = \Omega_+ \cup \Omega_-$, where
$\Omega_{\pm}$ are the disjoint open sets defined above, then Assumption 3.1 holds if and only if there is no non-trivial $v \in W_0^1(\Omega_\pm)$ that satisfies (6) in $\Omega_\pm$, which holds for all $k > 0$ outside a countable set whose only accumulation point is infinity.

In the case that Assumption 3.1 holds, and where $u$ is the unique solution to the above scattering problem and $u^t := u + u^i \in W_0^{1,\text{loc}}(\Omega)$, it proves convenient to extend $u^t$ by zero from $\Omega$ to $\mathbb{R}^n$ so that (as noted in §2.2) $u^t \in \tilde{H}^{1,\text{loc}}(\Omega) \subset H^{1,\text{loc}}(\mathbb{R}^n)$. We can correspondingly extend the definition of $u$ from $\Omega$ to $\mathbb{R}^n$, by setting $u := u^t - u^i \in H^{1,\text{loc}}(\mathbb{R}^n)$ (so that $u = -u^i$ on $\Gamma$ almost everywhere with respect to $n$-dimensional Lebesgue measure). We will assume these extensions hereafter, so that $u$ and $u^t$ are defined (almost everywhere) on $\mathbb{R}^n$ and $u, u^t \in H^{1,\text{loc}}(\mathbb{R}^n)$. Alternatively, one can require from the outset that $u \in H^{1,\text{loc}}(\mathbb{R}^n)$ and satisfies (6) in $\Omega$ and (7) and that $u^t = u + u^i \in \tilde{H}^{1,\text{loc}}(\Omega)$, in which case $u|_{\Omega}$ is the unique solution to the above scattering problem and $u = -u^i$ on $\Gamma$ (almost everywhere with respect to $n$-dimensional Lebesgue measure).

Introducing the orthogonal projection operator

$$ P: H^1(\mathbb{R}^n) \to \tilde{H}^1(\Omega)^\perp, \qquad (9) $$

we observe that $u^t \in \tilde{H}^{1,\text{loc}}(\Omega)$ if and only if $P(\sigma u^t) = 0$ for some, and hence every$^4$, $\sigma \in C_{0,\Gamma}^\infty$, in other words, if and only if

$$ P(\sigma u) = g, \qquad (10) $$

where

$$ g := -P(\sigma u^i). \qquad (11) $$

Thus (10) is an alternative formulation of the boundary condition that $u = -u^i$ on $\Gamma$, equivalent to the requirement that $u^t \in \tilde{H}^{1,\text{loc}}(\Omega)$.

To summarise, our scattering problem can be stated as follows: find $u \in H^{1,\text{loc}}(\mathbb{R}^n)$ satisfying (6) in $\Omega$, (7), and (10). We now reformulate this problem as an integral equation.

## 3.1 The integral equation on general compact sets

In what follows, $\mathcal{A}$ will denote the standard acoustic Newton potential operator, defined for compactly supported $\phi \in L_2(\mathbb{R}^n) = H^0(\mathbb{R}^n)$ by

$$ \mathcal{A}\phi(x) = \int_{\mathbb{R}^n} \Phi(x, y)\phi(y) \, dy, \quad x \in \mathbb{R}^n, \qquad (12) $$

<footnote>$^4$If $\sigma_1, \sigma_2 \in C_{0,\Gamma}^\infty$, then $\sigma_1 = \sigma_2$ on some open set $G \supset \Gamma$, so that, for $v \in H^{1,\text{loc}}(\mathbb{R}^n)$, $(\sigma_1 - \sigma_2)v \in H_{\mathbb{R}^n \setminus G}^{1,\text{loc}} \subset \tilde{H}^1(\Omega)$, so that $P((\sigma_1 - \sigma_2)v) = 0$.</footnote>
where $\Phi(x, y) := e^{ik|x-y|}/(4\pi|x-y|)$ ($n=3$), $\Phi(x, y) := \frac{i}{4}H_0^{(1)}(k|x-y|)$ ($n=2$), is the standard fundamental solution of the Helmholtz equation, and $H_0^{(1)}$ is the Hankel function of the first kind of order zero (e.g., [1, Eqn. (10.4.3)]). It is standard (see e.g. [33, Thm 3.1.2]) that, for $s \in \mathbb{R}$, in particular for $s=0$, $\mathcal{A}$ is continuous as a mapping

$$ \mathcal{A}: H_{\text{comp}}^{s-1}(\mathbb{R}^n) \to H^{s+1, \text{loc}}(\mathbb{R}^n), \quad (13) $$

where, for $s \in \mathbb{R}$, $H_{\text{comp}}^{s}(\mathbb{R}^n)$ is the space of compactly supported elements of $H^s(\mathbb{R}^n)$. Further (e.g. [33, Thm 3.1.4]),

$$ (\Delta + k^2)\mathcal{A}\phi = \mathcal{A}(\Delta + k^2)\phi = -\phi, \quad \phi \in H_{\text{comp}}^{-1}(\mathbb{R}^n). \quad (14) $$

Viewing $\mathcal{A}$ as an operator $\mathcal{A}: H_{\Gamma}^{-1} \to H^{1, \text{loc}}(\mathbb{R}^n) = W^{1, \text{loc}}(\mathbb{R}^n)$, we have that

$$ \mathcal{A}\phi(x) = \langle (\sigma\Phi(x, \cdot)), \bar{\phi} \rangle_{H^1(\mathbb{R}^n) \times H^{-1}(\mathbb{R}^n)}, \quad x \in \Omega, \quad (15) $$

where $\bar{\phi}$ denotes the complex conjugate of $\phi$ and $\sigma$ is any element of $C_{0,\Gamma}^{\infty}$ with $x \notin \text{supp } \sigma$. We define the operator $\mathcal{A}: H_{\Gamma}^{-1} \to \tilde{H}^1(\Omega)^{\perp} = (H_{\Gamma}^{-1})^*$ (the latter equality holding by [10, Cor. 3.4]) by

$$ \mathcal{A}\phi := P(\sigma\mathcal{A}\phi), \quad \phi \in H_{\Gamma}^{-1}, \quad (16) $$

with $\sigma \in C_{0,\Gamma}^{\infty}$ arbitrary. We also define the associated sesquilinear form $a(\cdot, \cdot)$ on $H_{\Gamma}^{-1} \times H_{\Gamma}^{-1}$ by

$$ a(\phi, \psi) := \langle A\phi, \psi \rangle_{H^1(\mathbb{R}^n) \times H^{-1}(\mathbb{R}^n)}, \quad \phi, \psi \in H_{\Gamma}^{-1}. \quad (17) $$

This form is compactly perturbed coercive, meaning that the operator $\mathcal{A}: H_{\Gamma}^{-1} \to \tilde{H}^1(\Omega)^{\perp}$ is a compact perturbation of a coercive operator (see, e.g., [11, §2.2] for detailed definitions and discussion). The following lemma is a generalisation of [13, Prop. 8.7, 8.8], and our proof is similar.

**Lemma 3.3.** Let $\Gamma \subset \mathbb{R}^n$ be compact. The sesquilinear form $a(\cdot, \cdot)$ is continuous and compactly perturbed coercive on $H_{\Gamma}^{-1} \times H_{\Gamma}^{-1}$, i.e., for some constants $C_a, \alpha > 0$, and some compact sesquilinear form $\tilde{a}(\cdot, \cdot)$,

$$ |a(\phi, \psi)| \le C_a \|\phi\|_{H_{\Gamma}^{-1}} \|\psi\|_{H_{\Gamma}^{-1}}, \quad |a(\phi, \phi) - \tilde{a}(\phi, \phi)| \ge \alpha \|\phi\|_{H_{\Gamma}^{-1}}^2, \quad \phi, \psi \in H_{\Gamma}^{-1}. \quad (18) $$

*Proof.* Continuity follows immediately from (13). Let us temporarily introduce the notations $\mathcal{A}_\kappa$ and $A_\kappa$ to denote the operators $\mathcal{A}$ and $A$ with $k$ replaced by some complex wavenumber $\kappa$. To prove that $a(\cdot, \cdot)$ is compactly perturbed coercive, we split the associated operator $\mathcal{A} = A_\kappa$ as $A_\kappa = A_i + (A_\kappa - A_i)$, where $A_i$ is the operator with wavenumber $\kappa = i$.
It is easy to check that for $\phi \in C_0^\infty(\mathbb{R}^n)$ the Fourier transform of $A_i\phi$ is given by $\hat{\phi}(\xi)/(|\xi|^2 + 1)$, which gives that

$$ \langle A_i\phi, \phi \rangle_{H^1(\mathbb{R}^n) \times H^{-1}(\mathbb{R}^n)} = \int_{\mathbb{R}^n} \frac{|\hat{\phi}(\xi)|^2}{|\xi|^2 + 1} d\xi = \|\phi\|_{H^{-1}(\mathbb{R}^n)}^2, \quad \phi \in H_{\Gamma}^{-1}, $$

so that $A_i : H_{\Gamma}^{-1} \to \tilde{H}^1(\Omega)^\perp$ is coercive with coercivity constant $\alpha = 1$. Further, for all $\sigma \in C_{0,\Gamma}^\infty$, $A_k - A_i = P\sigma(A_k - A_i)$ and, arguing as in [33, Rem. 3.1.3] and [13, Prop. 8.8], $(A_k - A_i) : H_{\Gamma}^{-1} \to H^{3,\text{loc}}(\mathbb{R}^n)$ is continuous, so that $\sigma(A_k - A_i) : H_{\Gamma}^{-1} \to H^1(\mathbb{R}^n)$ is compact, which implies that $A_k - A_i$ is compact. $\square$

Using Lemma 3.3 we can prove well-posedness of the scattering problem for arbitrary compact $\Gamma$, by reformulating it as a well-posed integral equation (IE), which we do in the following theorem. Our description of (20) as an IE will be justified in §3.2, when we discuss conditions under which the operator $A$ can be interpreted as an integral operator on $\Gamma$.

**Theorem 3.4.** Let $\Gamma \subset \mathbb{R}^n$ be compact, and suppose that Assumption 3.1 holds. Then $A : H_{\Gamma}^{-1} \to \tilde{H}^1(\Omega)^\perp$ is invertible. Further, for every $g \in \tilde{H}^1(\Omega)^\perp$ the problem defined by (6) in $\Omega$, (7), and (10) has a unique solution $u \in H^{1,\text{loc}}(\mathbb{R}^n)$ given by

$$ u = A\phi, \quad (19) $$

where $\phi \in H_{\Gamma}^{-1}$ is the unique solution of the IE

$$ A\phi = g, \quad (20) $$

which can be written equivalently in variational form as

$$ a(\phi, \psi) = \langle g, \psi \rangle_{H^1(\mathbb{R}^n) \times H^{-1}(\mathbb{R}^n)}, \quad \forall \psi \in H_{\Gamma}^{-1}. \quad (21) $$

If Assumption 3.1 does not hold then $A\phi = 0$ has a non-trivial solution $\phi \in H_{\Gamma}^{-1}$.

*Proof.* If $\phi \in H_{\Gamma}^{-1}$ satisfies (20) then $u$ given by (19) solves the scattering problem, by (14) and (16), and the fact that the acoustic Newton potential (12) satisfies (7). Thus if (20) has a solution then the scattering problem has a solution, and this solution is unique if Assumption 3.1 holds. By Lemma 3.3, $A$ is continuous and compactly perturbed coercive, and hence Fredholm of index zero by Lax-Milgram. Thus, to prove $A : H_{\Gamma}^{-1} \to \tilde{H}^1(\Omega)^\perp$ is invertible, so that (20) has a unique solution, it suffices to prove that $A$ is injective. For this, suppose that $\phi \in H_{\Gamma}^{-1}$ and $A\phi = 0$. Then $A\phi$ satisfies the homogeneous scattering problem, so, if Assumption 3.1 holds, we have $A\phi = 0$ in $\Omega$. But, for $\sigma \in C_{0,\Gamma}^\infty$, we also have $P(\sigma A\phi) = A\phi = 0$, so
$\sigma A\phi \in \tilde{H}^1(\Omega)$. Thus, by [10, eqn. (17)], $A\phi = \sigma A\phi = 0$ in $\Omega^c = \Gamma$, almost everywhere with respect to $n$-dimensional Lebesgue measure. Hence $A\phi = 0$, and by (14) we conclude that $\phi = 0$, proving injectivity, and hence invertibility of $A$.

If Assumption 3.1 does not hold then, by Remark 3.2, where $\Omega_-$ is as in that remark, there exists a non-zero $v \in \tilde{H}^1(\Omega_-) \subset \tilde{H}^1(\Omega)$ such that $\Delta v + k^2 v = 0$ in $\Omega_-$. By (14) we have $v = A\phi$, where $\phi := -(\Delta + k^2)v \in H_{\partial\Omega_-}^{-1} \subset H_{\Gamma}^{-1}$. Further, for $\sigma \in C_{0,\Gamma}^{\infty}$, $P(\sigma v) = 0$ since $\sigma v \in \tilde{H}^1(\Omega)$, so $A\phi = P(\sigma A\phi) = 0$, and $\phi \neq 0$ since $A\phi = v \neq 0$. $\square$

**Remark 3.5 (The role of the capacity of $\Gamma$).** We make the trivial observation that, if $H_{\Gamma}^{-1} = \{0\}$, then the only solution to (20) is $\phi = 0$, so that the scattered field $u = A\phi = 0$; i.e. the incident field does not interact with $\Gamma$. Further, $H_{\Gamma}^{-1} \neq \{0\}$ if and only if $\Gamma$ has positive $H^1$ capacity (see, e.g. [28, Thm 13.2.2], and for a collection of related results and generalisations, see [21, 10]). This holds if $\dim_H(\Gamma) > n - 2$ [21, Thm 2.12], and is equivalent to $\dim_H(\Gamma) > n - 2$ if $\Gamma$ is a $d$-set [21, Thm 2.17]. Moreover, by [10, Thm. 3.12], $H_{\Gamma}^{-1} = \{0\}$ if and only if $\tilde{H}^1(\Omega)^{\perp} = \{0\}$, so if $H_{\Gamma}^{-1} = \{0\}$ then the datum $g$ of the IE (20) is zero (recall the definition of $g$ in (9), (11)).

The following proposition concerns the support of the IE solution, and allows us to determine when the scattered fields and IE solutions for different scatterers $\Gamma$ coincide.

**Proposition 3.6.** Let $\Gamma \subset \mathbb{R}^n$ be compact with non-empty interior $\Gamma^\circ$ and let $\Omega := \Gamma^c$ be connected. Let $u = A\phi$ be the unique solution of the scattering problem for $\Gamma$, where $\phi \in H_{\Gamma}^{-1}$ is the unique solution of (20), with $g$ given by (11). Then $\phi \in H_{\partial\Gamma}^{-1}$. Suppose further that $\Gamma_\dagger$ is compact, with $\partial\Gamma \subset \Gamma_\dagger \subset \Gamma$, and that Assumption 3.1 is satisfied by $\Gamma_\dagger$, and let $u_\dagger = A\phi_\dagger$ be the unique solution of the scattering problem for $\Gamma_\dagger$, where $\phi_\dagger \in H_{\Gamma_\dagger}^{-1}$ is the unique solution of the IE for $\Gamma_\dagger$. Then $u_\dagger = u$ and $\phi_\dagger = \phi \in H_{\partial\Gamma}^{-1}$.

*Proof.* Since $A\phi + u^i = u + u^i = u^t \in \tilde{H}^{1,\text{loc}}(\Omega)$, it holds in $\Gamma^\circ$ that

$$0 = (\Delta + k^2)u^t = (\Delta + k^2)A\phi + (\Delta + k^2)u^i = -\phi, \quad (22)$$

by (14), and since $u^i$ satisfies (6) in a neighbourhood of $\Gamma$. Thus $\phi \in H_{\partial\Gamma}^{-1}$. Further, $u^t \in \tilde{H}^{1,\text{loc}}(\Omega)$ implies that $u^t = 0$, so $u = -u^i$, in $\Gamma^\circ$. Since $\Omega \subset \Omega_\dagger \subset \Omega \cup \Gamma^\circ$, where $\Omega_\dagger := \Gamma_\dagger^c$, and $u^i$ satisfies (6) in a neighbourhood of $\Gamma$, it follows that $u^t \in \tilde{H}^{1,\text{loc}}(\Omega_\dagger)$ and that $u$ satisfies (6) in $\Omega_\dagger$. Since $u$ also satisfies (7) it follows from Assumption 3.1 for $\Gamma_\dagger$ that $u_\dagger = u$. Since $u - u_\dagger = A(\phi - \phi_\dagger)$, it follows from (14) that $\phi = \phi_\dagger$. $\square$

**Remark 3.7.** Proposition 3.6 implies that if, for a given $k > 0$, Assumption 3.1 holds for a scatterer $\Gamma$ for which $\Omega := \Gamma^c$ is not connected, then the
scattered field and the IE solution $\phi$ for the scatterer $\Gamma$ coincide with those
for the scatterer $\Omega_+^c$, where $\Omega_+$ is the unbounded component of $\Gamma^c$, and $\phi$ is
supported in $\partial\Omega_+$.

**Remark 3.8 (Alternative IEs for the same scattering problem).**

Consider the case where $\Omega := \Gamma^c$ is connected and $\Gamma^\circ$ is non-empty (e.g., as in Figure 1(a) and (g)). By Proposition 3.6 and Remark 3.7, to solve the scattering problem for $\Gamma$ we can solve the IE on $\Gamma_\dagger$, for any compact $\Gamma_\dagger$ with $\partial\Gamma \subset \Gamma_\dagger \subset \Gamma$ provided $k^2$ is not a Dirichlet eigenvalue of $-\Delta$ in $\Omega_- := \Gamma \setminus \Gamma_\dagger$, in particular if $0 < k < k_0$ where $k_0^2$ is the smallest such eigenvalue. Recall from Theorem 3.4 that satisfying the IE on $\Gamma_\dagger$ is equivalent to requiring that $u^t = A\phi + u^i \in \tilde{H}^1(\Gamma_\dagger^c)$, i.e. to enforcing $A\phi = -u^i$ on $\Gamma_\dagger$.

Since $\phi \in H_{\partial\Gamma}^{-1}$, the choice $\Gamma_\dagger = \partial\Gamma$ is natural; with this choice the IE enforces $A\phi = -u^i$ on $\partial\Gamma$. When $\Omega$ is a Lipschitz open set this corresponds, as we discuss in Remark 3.17 below, to the standard single-layer-potential boundary IE (BIE) formulation. But it may be attractive to choose a larger $\Gamma_\dagger$, so that $A\phi = -u^i$ is enforced not just on $\partial\Gamma$ but also at points in $\Gamma^\circ$. This reduces the size of $\Omega_-$ and so increases $k_0$, and hence the interval $(0, k_0)$ in which the IE is uniquely solvable. (This is the rationale behind the CHIEF method and its variants for removing irregular frequencies of BIEs, e.g., [34], [37].) For the largest choice, $\Gamma_\dagger = \Gamma$, $A\phi = -u^i$ is enforced on the whole of $\Gamma$, and the IE is uniquely solvable for all $k > 0$, but at the cost in computation of discretising the whole of $\Gamma$ rather than $\partial\Gamma$ or some intermediate set.

We explore this further in §5, where we compare computations for the choices $\Gamma_\dagger = \Gamma$ and $\Gamma_\dagger = \partial\Gamma$ for the particular example of the Koch snowflake (Figure 1(g)) - see the discussion around Figure 7 below.

The variational formulation (21) will be the starting point for our Galerkin discretisation in §4. Having computed $\phi$ by solving a Galerkin discretisation of (21), we will evaluate $u(x)$ at points $x \in \Omega$ using the formulas (19)/(15). We will also compute the far-field pattern $u^\infty \in C^\infty(\mathbb{S}^{n-1})$, which satisfies (see, e.g., [7, Eqn. (2.23)], [29, p. 294])

$$u(x) = \frac{e^{ik|x|}}{|x|^{(n-1)/2}} (u^{\infty}(\hat{x}) + O(|x|^{-1})), \quad \text{as } |x| \to \infty,$$

uniformly in $\hat{x} := x/|x|$. Explicitly ([7, Eqn. (2.23)], [29, p. 294]),

$$u^{\infty}(\hat{x}) = \langle \sigma \Phi^{\infty}(\hat{x}, \cdot), \bar{\phi} \rangle_{H^1(\mathbb{R}^n) \times H^{-1}(\mathbb{R}^n)}, \quad \hat{x} \in \mathbb{S}^{n-1}, \quad (23)$$

where $\sigma$ is any element of $C_{0,\Gamma}^\infty$ and $\Phi^\infty(\cdot, y)$ is the far-field pattern of $\Phi(\cdot, y)$, for $y \in \mathbb{R}^n$, viz.

$$\Phi^{\infty}(\hat{x}, y) := \frac{ik^{(n-3)/2}}{2(2\pi i)^{(n-1)/2}} \exp(-ik\hat{x} \cdot y), \quad \hat{x} \in \mathbb{S}^{n-1}, y \in \mathbb{R}^n. \quad (24)$$
## 3.2 The integral equation on $d$-sets in trace spaces

So far, our analysis has been for general compact scatterers $\Gamma \subset \mathbb{R}^n$. We now assume additionally that $\Gamma$ is a $d$-set (in the sense of (2)), and that $n - 2 < d \le n$, so we have a non-trivial scattered field (see Remark 3.5). In this case one can view the operator $A$ as an integral operator with respect to Hausdorff measure $\mathcal{H}^d$, by reinterpreting $A$ as a map between certain “trace spaces” on $\Gamma$. This will allow us to relate our IE (20) to previously studied IE formulations in certain special cases (see Remark 3.17), and will pave the way for the discretization we consider in §4. We begin by briefly recalling the definition of trace spaces on $d$-sets, and the relationship between them and function spaces on $\mathbb{R}^n$. For a more detailed explanation see [5, §2.4].

Let $\Gamma \subset \mathbb{R}^n$ be a $d$-set for some $n - 2 < d \le n$.<fnref n="5" /> We denote by $L_2(\Gamma)$ the Hilbert space of functions on $\Gamma$ that are measurable and square integrable with respect to $\mathcal{H}^d|_{\Gamma}$, normed by $\|f\|_{L_2(\Gamma)} := (\int_{\Gamma} |f(x)|^2 d\mathcal{H}^d(x))^{1/2}$, and by $L_\infty(\Gamma)$ the Banach space of functions on $\Gamma$ that are measurable and essentially bounded with respect to $\mathcal{H}^d|_{\Gamma}$, normed by $\|f\|_{L_\infty(\Gamma)} := \ess\sup_{x \in \Gamma} |f(x)|$.

Let $\text{tr}_\Gamma : C_0^\infty(\mathbb{R}^n) \to L_2(\Gamma)$ be the trace (or restriction) operator, with dense range, defined by $\text{tr}_\Gamma\varphi = \varphi|_{\Gamma} \in L_2(\Gamma)$, for $\varphi \in C_0^\infty(\mathbb{R}^n)$. For $s > \frac{n-d}{2}$, this extends to a continuous linear operator

$$ \text{tr}_{\Gamma}: H^s(\mathbb{R}^n) \to L_2(\Gamma), \qquad (25) $$

also with dense range (see [36, Thm 18.6] and [5, §2.4]). Setting

$$ t := s - \frac{n-d}{2} > 0, \qquad (26) $$

we define the trace space $\mathcal{H}^t(\Gamma) := \text{tr}_\Gamma(H^s(\mathbb{R}^n)) \subset L_2(\Gamma)$, and equip it with the quotient norm

$$ \|f\|_{\mathcal{H}^t(\Gamma)} := \inf_{\substack{\varphi \in H^s(\mathbb{R}^n) \\ \text{tr}_\Gamma \varphi = f}} \|\varphi\|_{H^s(\mathbb{R}^n)}. $$

This makes $\mathcal{H}^t(\Gamma)$ a Hilbert space unitarily isomorphic to the quotient space $H^s(\mathbb{R}^n)/\ker(\text{tr}_\Gamma)$. For $t > 0$ we denote by $\mathcal{H}^{-t}(\Gamma)$ the dual space $(\mathcal{H}^t(\Gamma))^*$. Identifying $L_2(\Gamma)$ with its dual in the standard way, and with $\mathcal{H}^0(\Gamma) := L_2(\Gamma)$, we have that $\mathcal{H}^{t'}(\Gamma)$ is continuously embedded in $\mathcal{H}^t(\Gamma)$ with dense image for any $t, t' \in \mathbb{R}$ with $t' > t$, and if $g \in \mathcal{H}^t(\Gamma)$ for some $t \ge 0$ and $f \in L_2(\Gamma)$ then

$$ \langle f, g \rangle_{\mathcal{H}^{-t}(\Gamma) \times \mathcal{H}^t(\Gamma)} = (f, g)_{L_2(\Gamma)}. \qquad (27) $$

<footnote n="5">⁵While our focus here is on $n = 2, 3$ and $n - 2 < d \le n$, we note that the definitions and results from the current paragraph onwards, up to and including Lemma 4.1, all extend to general $n \in \mathbb{N}$ and $0 < d \le n$. For details see, e.g., [5, §2.4].</footnote>
Assuming (26), $\text{tr}_\Gamma : H^s(\mathbb{R}^n) \to \mathbb{H}^t(\Gamma)$ is a continuous linear surjection that has unit norm and is a unitary isomorphism from $\text{ker}(\text{tr}_\Gamma)^\perp$ to $\mathbb{H}^t(\Gamma)$. Furthermore, $\tilde{H}^s(\Gamma^c) \subset \text{ker}(\text{tr}_\Gamma)$. Accordingly, again assuming (26), its adjoint

$$ \text{tr}_{\Gamma}^{*} : \mathbb{H}^{-t}(\Gamma) \to H^{-s}(\mathbb{R}^{n}), \qquad (28) $$

is a continuous linear isometry with range contained in $H_{\Gamma}^{-s} = (\tilde{H}^s(\Gamma^c)^{\perp})^*$, satisfying

$$ \langle \varphi, \text{tr}_{\Gamma}^{*} f \rangle_{H^s(\mathbb{R}^n) \times H^{-s}(\mathbb{R}^n)} = \langle \text{tr}_{\Gamma} \varphi, f \rangle_{\mathbb{H}^t(\Gamma) \times \mathbb{H}^{-t}(\Gamma)}, \quad f \in \mathbb{H}^{-t}(\Gamma), \varphi \in H^s(\mathbb{R}^n). \qquad (29) $$

In particular, when $f \in \mathbb{L}_2(\Gamma)$ we have that

$$ \langle \varphi, \text{tr}_{\Gamma}^{*} f \rangle_{H^s(\mathbb{R}^n) \times H^{-s}(\mathbb{R}^n)} = (\text{tr}_{\Gamma} \varphi, f)_{\mathbb{L}_2(\Gamma)}. \qquad (30) $$

The following theorem is a slight modification of results presented in [3] and [5].

**Theorem 3.9.** Let $\Gamma \subset \mathbb{R}^n$ be a $d$-set for some $n-2 < d \le n$. If $\frac{n-d}{2} < s < \frac{n-d}{2} + 1$, so that $t := s - \frac{n-d}{2} \in (0,1)$, then for $\text{tr}_\Gamma : H^s(\mathbb{R}^n) \to \mathbb{H}^t(\Gamma)$ it holds that $\text{ker}(\text{tr}_\Gamma) = \tilde{H}^s(\Gamma^c)$.

Hence $\text{tr}_\Gamma|_{\tilde{H}^s(\Gamma^c)^\perp} : \tilde{H}^s(\Gamma^c)^\perp \to \mathbb{H}^t(\Gamma)$ is a unitary isomorphism, the range of $\text{tr}_\Gamma^*$ is equal to $H_\Gamma^{-s}$, the map $\text{tr}_\Gamma^* : \mathbb{H}^{-t}(\Gamma) \to H_\Gamma^{-s}$ is a unitary isomorphism, and $\text{tr}_\Gamma^*(\mathbb{L}_2(\Gamma))$ is dense in $H_\Gamma^{-s}$.

If $d=n$ (in which case $t=s$) then the above statements hold also for the limiting case $s=t=0$.

*Proof.* For the statements for $\frac{n-d}{2} < s < \frac{n-d}{2} + 1$ see [3, Prop. 6.7, Thm 6.13] and [5, Thm 2.7]). For the limiting case mentioned, we note that when $d=n$ we have $\mathbb{L}_2(\Gamma) = L_2(\Gamma)$, and the trace map $\text{tr}_\Gamma : H^s(\mathbb{R}^n) \to \mathbb{L}_2(\Gamma)$ is continuous for all $s \ge 0$ and is given simply by $\text{tr}_\Gamma \phi = \phi|_\Gamma$, so that $\text{ker}(\text{tr}_\Gamma)$ is the set of functions in $H^s(\mathbb{R}^d)$ that vanish almost everywhere on $\Gamma$ with respect to Lebesgue measure. For $s=0$, that $\text{ker}(\text{tr}_\Gamma) = \tilde{H}^0(\Gamma^c)$ is then a consequence of [10, Lem. 3.16]. $\square$

**Remark 3.10 (Connection to known cases: I).** The spaces $\mathbb{H}^t(\Gamma)$ introduced above can be related to well-known trace spaces in special cases. For instance:

(a) If $\Gamma$ is the closure of a bounded Lipschitz open set $D$ (e.g. Figure 1(a)), then $\Gamma$ is a $d$-set with $d=n$, so that $t=s$, and $\mathbb{H}^t(\Gamma)$ coincides with the restriction space $H^t(D) = \{U|_\Gamma : U \in H^t(\mathbb{R}^n)\}$ for $t \ge 0$ (this follows from the fact that $\text{tr}_\Gamma$ coincides with the restriction operator $|_\Gamma$ in this case, as noted in the proof of Theorem 3.9 above). Hence, for $t < 0$, $\mathbb{H}^t(\Gamma)$ is unitarily isomorphic to the space $\tilde{H}^t(D) = (H^{-t}(D))^*$ (see, e.g., [29, Thm 3.30]).
(b) If $\Gamma$ is the boundary of a bounded Lipschitz open set (e.g. Figure 1(b)), then $\Gamma$ is a $d$-set with $d = n-1$ and $\mathbb{H}^t(\Gamma)$ coincides, for $|t| \le 1/2$, with the boundary Sobolev space $H^t(\Gamma)$, as defined, e.g., in [29, pp. 98-99]; see [3, Rem. 6.5] for details.

(c) A further example of a family of $d$-sets with $d = n-1$ is provided by the “multi-screens” defined in Definition 2.3 of [13]; these are finite unions of Lipschitz subsets of the boundaries of bounded Lipschitz open sets, a specific example being given in Figure 1(d). In [13], trace spaces on multi-screens are defined as quotient spaces. Specifically, in the parlance of multi-screen theory, the space $H^{1/2}([\Gamma]) := H^1(\mathbb{R}^n)/\tilde{H}^1(\Gamma^c)$ is referred to as the “Dirichlet single-trace space” (see [13, Defn 6.1]). It is unitarily isomorphic to the space $\mathbb{H}^{1/2}(\Gamma)$. The dual space $\tilde{H}^{-1/2}([\Gamma]) := (H^{1/2}([\Gamma]))^*$ is referred to as the “Neumann jump space” (see [13, Defn 6.4]). It is unitarily isomorphic to $\mathbb{H}^{-1/2}(\Gamma) = (\mathbb{H}^{1/2}(\Gamma))^*$, which, by Theorem 3.9, is unitarily isomorphic to $H_{\Gamma}^{-1}$.

In [5] we studied scattering by planar screens $\Gamma \subset \Gamma_\infty := \mathbb{R}^{n-1} \times \{0\}$ (see Figure 1(c) and (e)). There we defined trace spaces on $\Gamma$ by a two-step process, first taking a trace onto the hyperplane $\Gamma_\infty$, then applying the above trace results in $\mathbb{R}^{n-1}$. The next result, which is stated without proof since it follows trivially from standard trace mapping properties (see e.g. [29, Lem. 3.35]), allows us to avoid this complication, and instead treat a planar screen as any other compact set.

**Lemma 3.11 (Planar screens).** Suppose that $\Gamma = \tilde{\Gamma} \times \{0\} \subset \Gamma_\infty := \mathbb{R}^{n-1} \times \{0\}$ where $\tilde{\Gamma} \subset \mathbb{R}^{n-1}$ is a $d$-set in $\mathbb{R}^{n-1}$ for some $n-2 < d \le n-1$. Then $\Gamma$ is a $d$-set in $\mathbb{R}^n$. Further, noting from (26) that $t = s - (n-d)/2 = (s - 1/2) - (n-1-d)/2$, for $s \ge (n-d)/2$ we have one continuous trace operator $\text{tr}_{\tilde{\Gamma}}: H^{s-1/2}(\mathbb{R}^{n-1}) \to L_2(\tilde{\Gamma})$ with range $\mathbb{H}^t(\tilde{\Gamma})$, and another continuous trace operator $\text{tr}_{\Gamma}: H^s(\mathbb{R}^n) \to L_2(\Gamma)$ with range $\mathbb{H}^t(\Gamma)$. The spaces $\mathbb{H}^t(\Gamma)$ and $\mathbb{H}^t(\tilde{\Gamma})$ coincide up to identification by the map $f \mapsto f(\cdot, 0)$, and the traces satisfy $\text{tr}_{\Gamma} = \text{tr}_{\tilde{\Gamma}} \circ \gamma$ where $\gamma: H^s(\mathbb{R}^n) \to H^{s-1/2}(\mathbb{R}^{n-1})$ is the standard (surjective) trace operator.

Of central importance is the case $s = 1$, as our scattering problem is posed in $H^{1,\text{loc}}(\mathbb{R}^n)$. Hence we give the value of $t$ given by (26) for $s = 1$ its own notation, defining (for $n-2 < d \le n$)

$$t_d := 1 - \frac{n-d}{2} \in (0, 1]. \qquad (31)$$

The case $s = 1$ ($t = t_d$) is covered by Theorem 3.9 for $n-2 < d < n$, but not for $d = n$. It is not known to us whether $\text{ker}(\text{tr}_\Gamma) = \tilde{H}^1(\Gamma^c)$ for general $n$-sets $\Gamma$, although we know it to hold in many cases (see Remark 3.14).
<figure><img src="image_4.png" /><figcaption>Figure 3: Schema of relevant function spaces and operators for $s = 1$, $t = t_d := 1 - \frac{n-d}{2} \in (0, 1]$. The operators $\text{tr}_\Gamma$ and $\text{tr}_\Gamma^*$ are isometries, indeed unitary isomorphisms if Assumption 3.12 holds.</figcaption></figure>

For convenience we introduce this as an assumption, with which assumption Corollary 3.13 below is an immediate consequence of Theorem 3.9.

**Assumption 3.12.** $\Gamma \subset \mathbb{R}^n$ is a $d$-set with either (i) $n - 2 < d < n$ or (ii) $d = n$ and $\text{ker}(\text{tr}_\Gamma) = \tilde{H}^1(\Gamma^c)$.

**Corollary 3.13.** Suppose that Assumption 3.12 holds. Then, with $t_d$ defined by (31), $\text{tr}_\Gamma : \tilde{H}^1(\Gamma^c)^\perp \to \mathbb{H}^{t_d}(\Gamma)$ and $\text{tr}_\Gamma^* : \mathbb{H}^{-t_d}(\Gamma) \to H_\Gamma^{-1}$ are unitary isomorphisms and $\text{tr}_\Gamma^*(\mathbb{L}_2(\Gamma))$ is dense in $H_\Gamma^{-1}$.

**Remark 3.14.** In the case $d = n$, a sufficient condition for $\text{ker}(\text{tr}_\Gamma) = \tilde{H}^1(\Gamma^c)$ is that $\tilde{H}^1(\Gamma^c) = H_{\Gamma^c}^{1,0}$, since $\tilde{H}^1(\Gamma^c) \subset \text{ker}(\text{tr}_\Gamma) \subset H_{\Gamma^c}^{1,0}$ (see, e.g., [10, Eqn (17)]). Hence Assumption 3.12 holds in the following cases:

(i) $\Gamma^c$ is $C^0$ [29, Thm. 3.29], so in particular if $\Gamma$ is the closure of a Lipschitz open set; more generally, if $\Gamma^c$ is $C^0$ except at points in a closed, countable subset $P$ of $\partial(\Gamma^c)$ with at most finitely many limit points [10, Thm. 3.24]; see, e.g., the examples in [10, Fig. 4];

(iii) $\Gamma$ is an $n$-set OSC-IFS attractor (see [4]), e.g. the Koch snowflake in Figure 1(g);

(iv) $\Gamma$ is the closure of one of the classical snowflake domains in [3, §5.1].

Returning to our integral equation, we can now “lift” the potential $\mathcal{A}: H_\Gamma^{-1} \to H^{1,\text{loc}}(\mathbb{R}^n)$ and the operator $A: H_\Gamma^{-1} \to \tilde{H}^1(\Omega)^\perp$ to continuous maps on the trace spaces via the compositions

$$ \begin{gathered} \mathcal{A}\text{tr}_\Gamma^* : \mathbb{H}^{-t_d}(\Gamma) \to H^{1,\text{loc}}(\mathbb{R}^n) \quad \text{and} \\ \mathcal{A} := \text{tr}_\Gamma A \text{tr}_\Gamma^* : \mathbb{H}^{-t_d}(\Gamma) \to \mathbb{H}^{t_d}(\Gamma). \end{gathered} \tag{32} $$
The next lemma collects some basic results about these compositions. A schematic showing the relationships between the main function spaces and operators involved is given in Figure 3.

**Lemma 3.15.** Let $\Gamma$ be a compact $d$-set with $n - 2 < d \le n$, and define $\mathcal{A}$ as in (32). Then:

$$ (i) \quad \text{For arbitrary } \sigma \in C_{0,\Gamma}^{\infty}, \qquad \mathcal{A}\Psi = \operatorname{tr}_{\Gamma}(\sigma \operatorname{Atr}_{\Gamma}^{*}\Psi), \quad \Psi \in \mathbb{H}^{-t_d}(\Gamma); \qquad (33) $$

(ii) The sesquilinear form on $\mathbb{H}^{-t_d}(\Gamma) \times \mathbb{H}^{-t_d}(\Gamma)$ associated with the operator $\mathcal{A}$ satisfies

$$ \langle \mathcal{A}\Psi, \tilde{\Psi} \rangle_{\mathbb{H}^{t_d}(\Gamma) \times \mathbb{H}^{-t_d}(\Gamma)} = a(\operatorname{tr}_{\Gamma}^{*}\Psi, \operatorname{tr}_{\Gamma}^{*}\tilde{\Psi}), \quad \Psi, \tilde{\Psi} \in \mathbb{H}^{-t_d}(\Gamma), \qquad (34) $$

with $a(\cdot, \cdot)$ as in (17), and is continuous and compactly perturbed co-ercive;

(iii) If Assumption 3.12 holds, then problem (20)/(21) can be equivalently stated as: given $g \in \tilde{H}^1(\Omega)^{\perp}$, find $\Psi \in \mathbb{H}^{-t_d}(\Gamma)$ such that

$$ \mathcal{A}\Psi = \operatorname{tr}_{\Gamma} g, \qquad (35) $$

or, equivalently,

$$ \langle \mathcal{A}\Psi, \tilde{\Psi} \rangle_{\mathbb{H}^{t_d}(\Gamma) \times \mathbb{H}^{-t_d}(\Gamma)} = \langle \operatorname{tr}_{\Gamma} g, \tilde{\Psi} \rangle_{\mathbb{H}^{t_d}(\Gamma) \times \mathbb{H}^{-t_d}(\Gamma)}, \quad \tilde{\Psi} \in \mathbb{H}^{-t_d}(\Gamma), \qquad (36) $$

with the solutions of (20)/(21) and (35)/(36) related by $\phi = \operatorname{tr}_{\Gamma}^{*}\Psi$. If Assumption 3.1 also holds, then $A : H_{\Gamma}^{-1} \to \tilde{H}^1(\Omega)^{\perp}$ and $\mathcal{A} : \mathbb{H}^{-t_d}(\Gamma) \to \mathbb{H}^{t_d}(\Gamma)$ are both invertible.

*Proof.* (i) Equation (33) follows from (9), (16), and the fact that $\tilde{H}^1(\Omega) \subset \ker(\operatorname{tr}_{\Gamma})$, which implies that $\operatorname{tr}_{\Gamma} P\phi = \operatorname{tr}_{\Gamma} \phi$ for $\phi \in H^1(\mathbb{R}^n)$.

(ii) Equation (34) follows from (17) and (29), and the rest of (ii) follows from (34), Lemma 3.3, and that $\operatorname{tr}_{\Gamma}^{*} : \mathbb{H}^{-t_d}(\Gamma) \to H_{\Gamma}^{-1}$ is an isometry.

(iii) The first statement follows from Corollary 3.13, and the second from Theorem 3.4. $\square$

Crucial for the practical implementation of the Hausdorff IE method described in §4 is the fact that both $\operatorname{Atr}_{\Gamma}^{*}$ and $\mathcal{A}$ have integral representations with respect to Hausdorff measure.

**Theorem 3.16.** Let $\Gamma$ be a compact $d$-set with $n - 2 < d \le n$. Then:

$$ (i) \quad \text{For } \Psi \in \mathbb{L}_2(\Gamma), \qquad \operatorname{Atr}_{\Gamma}^{*}\Psi(x) = \int_{\Gamma} \Phi(x, y)\Psi(y) \, d\mathcal{H}^d(y), \quad x \in \Omega; \qquad (37) $$
(ii) For $\Psi \in L_{\infty}(\Gamma)$, the right-hand side of (37) is well-defined (as a Lebesgue integral with respect to $H^d$ measure) for all $x \in \mathbb{R}^n$, and is a continuous function on $\mathbb{R}^n$. Further, (37) holds for almost all $x \in \mathbb{R}^n$ with respect to n-dimensional Lebesgue measure, so that $\mathrm{Atr}_{\Gamma}^*\Psi \in C(\mathbb{R}^n)$.

$$ (iii) \text{ For } \Psi \in L_{\infty}(\Gamma), \\ \qquad \mathbb{A}\Psi(x) = \int_{\Gamma} \Phi(x, y)\Psi(y)d\mathcal{H}^d(y), \quad \text{for } \mathcal{H}^d\text{-a.e. } x \in \Gamma. \quad (38) $$

*Proof.* (i) This is an immediate consequence of (15) and (30).

(ii) For $\Psi \in L_{\infty}(\Gamma)$, that the right-hand-side of (37) is well-defined for all $x \in \mathbb{R}^n$ and defines a continuous function on $\mathbb{R}^n$ follows as in the proof of [5, Prop. 4.5], using the estimates for convolution integrals with respect to $H^d$ measure on $d$-sets in [5, Rem. 2.2] (cf. [26, Lemma 2.18]). For $d < n$, in which case $m(\partial\Gamma) = H^n(\Gamma) = 0$, we have from (i) that (37) holds for almost all $x \in \mathbb{R}^n$. For $d = n$, when $\Psi \in L_{\infty}(\Gamma) \subset L_2(\mathbb{R}^n)$, that (37) holds for almost all $x \in \mathbb{R}^n$ is just a special case of (12).

(iii) Suppose $\Psi \in L_{\infty}(\Gamma)$. For $x \in \mathbb{R}^n$ let $G(x)$ denote the right-hand side of (37). By part (ii) and (13), $\sigma G \in H^1(\mathbb{R}^n)$ and is continuous, for $\sigma \in C_{0,\Gamma}^{\infty}$, so that $\mathrm{tr}_{\Gamma}(\sigma G) = (\sigma G)|_{\Gamma} = G|_{\Gamma}$. Since $\mathrm{Atr}_{\Gamma}^*\Psi(x) = G(x)$, for a.e. $x \in \mathbb{R}^n$ with respect to n-dimensional Lebesgue measure by part (ii), it follows by (33) that $\mathbb{A}\Psi = G|_{\Gamma}$ in $L_2(\Gamma)$, giving the claimed result. $\square$

Our next remark builds on the characterisations of the spaces $H^l(\Gamma)$ in Remark 3.10.

**Remark 3.17 (Connection to known cases: II).** The trace space formulation (35) of our IE is familiar in a number of cases, in each of which $\Gamma$ is a $d$-set by Remark 3.10 or Lemma 3.11.

(a) If $\Gamma$ is the boundary of a bounded Lipschitz open set (see, e.g., Figure 1(b) and Remark 3.10(b)), then $d = n-1$, $t_d = 1/2$, and $H^d$ coincides with surface measure on $\Gamma$ [15, Theorem 3.8]. Thus the expression (37) for $\mathrm{Atr}_{\Gamma}^*\Psi$ coincides with the definition (e.g., [7, Eqn. (2.19)]) of the standard single-layer potential with density $\Psi$. The representation (38) for $\mathbb{A}$ coincides with the definition [7, Eqn. (2.32)] of the single-layer boundary integral operator $S: H^{-1/2}(\Gamma) \to H^{1/2}(\Gamma)$, viz. $S\phi(x) = \int_{\Gamma} \Phi(x, y)\phi(x)ds(y)$, and (33) coincides with [33, Eqn (3.6), Def. 3.15]. The IE (35) coincides with [7, Eqn. (2.63)].

(b) In the case where $\Gamma$ is a multi-screen (see, e.g., Figure 1(d) and Remark 3.10(c)), $d = n-1$, $t_d = 1/2$, the representation (37) for $\mathrm{Atr}_{\Gamma}^*$ coincides with the definition of the single-layer potential given in [13, Eqn. (8.2)] and the representation (38) for $\mathbb{A}$ coincides with the first
boundary integral operator from Proposition 8.8 in [13] (with the same explicit surface integral form as in part (a) above). The mapping properties and coercivity up to a compact perturbation derived in Lemma 3.15(ii) generalize the first inequality of [13, Prop. 8.8].

(c) If the $d$-set $\Gamma$ is a compact subset of $\Gamma_\infty := \mathbb{R}^{n-1} \times 0$ (see Lemma 3.11), so $\Gamma$ is a planar screen (examples are Figure 1(c) and (e)), the expression (37) for $\mathrm{Atr}_\Gamma^* \Psi$ coincides with [5, Eqn. (43)] and the representation (38) for $\mathrm{A}$ coincides with [5, Eqn. (49)].

The definition and mapping properties of $\mathrm{tr}_\Gamma$ and $\mathrm{tr}_\Gamma^*$, noted in (25) and (28), combined with the representation (33), enable us to extend the domain of $\mathrm{A}$ to $\mathbb{H}^{-t}(\Gamma)$ for $t_d < t < 2t_d$ or restrict it to $\mathbb{H}^{-t}(\Gamma)$ for $0 < t < t_d$ as stated in the following result (cf. [5, Prop. 4.7]).

**Proposition 3.18.** Let $\Gamma$ be a compact $d$-set with $n-2 < d \le n$, and let $|t| < t_d$. Then $\mathrm{A} : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ and is continuous. When $d=n$ this holds also for $t = \pm t_d = \pm 1$.

*Proof.* This follows from (33) and the mapping properties (13) of $\mathcal{A}$, recalling from (25) and (28) that, for $s > (n-d)/2$, the mappings $\mathrm{tr}_\Gamma : H^s(\mathbb{R}^n) \to \mathbb{H}^t(\Gamma)$ and $\mathrm{tr}_\Gamma^* : \mathbb{H}^{-t}(\Gamma) \to H_\Gamma^{-s}$ are continuous, where $t = s - (n-d)/2 > 0$, and that, when $d=n$ and $t=s$, these mappings are continuous also for $s=0$ (see the proof of Theorem 3.9). $\square$

Under certain assumptions, $\mathrm{A} : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ is invertible for a range of $t$ around 0.

**Proposition 3.19.** Let $\Gamma$ be an OSC-IFS attractor with dimension $d := \dim_{\mathbb{H}}(\Gamma)$ such that either (a) $\Gamma$ is disjoint with $n-2 < d < n$, (b) $d=n$, or (c) $d=n-1$ and $\Gamma \subset \Gamma_\infty := \mathbb{R}^{n-1} \times \{0\}$. If Assumption 3.1 holds, there exists $0 < \epsilon \le t_d$ such that $\mathrm{A} : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ is invertible for $|t| < \epsilon$.

*Proof.* The claimed invertibility of $\mathrm{A} : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ for a range of $t$ in a neighbourhood of $t=0$ follows by applying a result on interpolation of invertibility of operators ([30, Prop. 4.7], which quotes [35]), recalling that (i) $\mathrm{A} : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ is bounded for $|t| < t_d$ (Proposition 3.18); (ii) $\mathrm{A} : \mathbb{H}^{-t_d}(\Gamma) \to \mathbb{H}^{t_d}(\Gamma)$ is invertible, as noted below (36); and (iii) in the case that $\Gamma$ is disjoint and $d < n$, $\{\mathbb{H}^t(\Gamma)\}_{|t|<1}$ is an interpolation scale [5, Cor. 3.3]; (iv) in the case that $d=n$, $\{\mathbb{H}^t(\Gamma)\}_{t\ge0}$ and $\{\mathbb{H}^t(\Gamma)\}_{t\le0}$ are interpolation scales [4]; (v) in the case $d=n-1$ and $\Gamma \subset \Gamma_\infty := \mathbb{R}^{n-1} \times \{0\}$ the statement of (iv) is also true, by (iv) applied in $\mathbb{R}^{n-1}$ and Lemma 3.11. $\square$

Invertibility of $\mathrm{A} : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ for some $t > 0$ implies a regularity result about the solution $\phi$ of the IE (20), provided that the datum $g$ is sufficiently smooth.
**Remark 3.20 (Solution regularity in the $H_Γ^s$ scale).** If there exists $0 < t < t_d$ such that $A : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ is invertible and $\text{tr}_\Gamma g \in \mathbb{H}^{t+t_d}(\Gamma)$, and if Assumption 3.12 holds, then, by (32) and the mapping properties of $\text{tr}_\Gamma^*$ recalled in Theorem 3.9, the solution $\phi = \text{tr}_\Gamma^* \Psi$ of the IE (20) satisfies

$$ \phi \in H_{\Gamma}^{-1+t}, \quad \text{with} \quad \|\phi\|_{H_{\Gamma}^{-1+t}} \le C \|\text{tr}_{\Gamma} g\|_{\mathbb{H}^{t+t_d}(\Gamma)}, \qquad (39) $$

for some constant $C > 0$ independent of $\phi$ and $g$. In the case of scattering of an incident wave $u^i$, in which $g$ is given by (11), we have that $\text{tr}_\Gamma g = -u^i|_\Gamma \in \mathbb{H}^{t+t_d}(\Gamma)$ for all $0 < t < t_d$, since $u^i$ is $C^\infty$ in a neighbourhood of $\Gamma$. Hence, in this case, if the conditions of Proposition 3.19 hold, then (39) holds for $0 < t < \epsilon$, where $\epsilon$ is as in Proposition 3.19.

Given $\Gamma$, determining the largest value of $t$ for which $A : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ is invertible is an open problem. However, so that we have a theoretical prediction against which to compare our numerical results in §5, we consider the following hypothesis.

**Hypothesis 3.21.** $A : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ is invertible for $0 \le t < t_d'$, where $t_d' := 1 - (n-d')/2$, with $d' := \dim_H(\partial\Omega_+)$ and $\Omega_+$ the unbounded component of $\Omega = \Gamma^c$.

To give some context for Hypothesis 3.21, we note that, to match Proposition 3.18, a naive hypothesis might be that $A : \mathbb{H}^{t-t_d}(\Gamma) \to \mathbb{H}^{t+t_d}(\Gamma)$ is invertible for all $0 \le t < t_d$ (cf. [5, Conj. 4.8] in the planar screen case). However, such a hypothesis fails in cases where $d' < d$, such as Figures 1(a) (where $d' = 1 < 2 = d$) and 1(g) (where $d' = \log(4)/\log(3) < 2 = d$). Indeed, if this naive hypothesis were to hold, then, for a scattering problem with $u^i|_\Gamma \ne 0$, arguing as in Remark 3.20 it would follow that $0 \ne \phi \in H_\Gamma^s$ for every $s < -(n-d)/2$. Furthermore, by Remark 3.7 we would have that $\text{supp}\phi \subset \partial\Omega_+$, from which it would follow that $0 \ne \phi \in H_{\partial\Omega_+}^s$ for every $s < -(n-d)/2$. But if $d' < d$ this is impossible because $H_{\partial\Omega_+}^s = \{0\}$ for $-(n-d')/2 < s < -(n-d)/2$, in fact for $-(n-d')/2 \le s < -(n-d)/2$ if $\partial\Omega_+$ is a $d'$-set [21, Thms 2.12 & 2.17]. Therefore, Hypothesis 3.21 is the strongest hypothesis that is consistent with Remark 3.7. In §5 we report numerical results which suggest that Hypothesis 3.21 may hold in certain cases, but not in general.

# 4 The Hausdorff-measure IEM

We now define and analyse our Hausdorff-measure Galerkin IEM. To begin with, let us assume that $\Gamma$ is a compact $d$-set for some $n-2 < d \le n$. Given $N \in \mathbb{N}$ let $\{T_j\}_{j=1}^N$ be a mesh of $\Gamma$, a collection of $\mathcal{H}^d$-measurable subsets of
$\Gamma$ (the elements) such that

$$ \Gamma = \bigcup_{j=1}^{N} T_j, \quad \mathcal{H}^d(T_j) > 0 \text{ for } j = 1, \dots, N, \quad \text{and} \quad \mathcal{H}^d(T_j \cap T_{j'}) = 0 \text{ for } j \neq j', $$

and set $h := \max_{j=1,\dots,N} \text{diam}(T_j)$. Define the $N$-dimensional space of piecewise constants

$$ \mathbb{V}_N := \{f \in \mathbb{L}_2(\Gamma) : f|_{T_j} = c_j \text{ for some } c_j \in \mathbb{C}, j = 1, \dots, N\} \subset \mathbb{L}_2(\Gamma) \quad (40) $$

and set

$$ V_N := \text{tr}_{\Gamma}^*(\mathbb{V}_N) \subset H_{\Gamma}^{-1}. \quad (41) $$

Under appropriate assumptions, the spaces $V_N$ are dense in $H_{\Gamma}^{-1}$ as $N \to \infty$.

**Lemma 4.1.** Suppose that Assumption 3.12 holds, and that $h \to 0$ as $N \to \infty$. Then

$$ \inf_{\psi_N \in V_N} \|\psi - \psi_N\|_{H^{-1}(\mathbb{R}^n)} \to 0 \quad \text{as} \quad N \to \infty, \quad \text{for all} \quad \psi \in H_{\Gamma}^{-1}. \quad (42) $$

*Proof.* Suppose that $h \to 0$ as $N \to \infty$. It is easy to see (see the proof of [5, Thm 5.1]) that $\inf_{f_N \in \mathbb{V}_N} \|f - f_N\|_{\mathbb{L}_2(\Gamma)} \to 0$ as $N \to \infty$ for every $f \in \mathbb{L}_2(\Gamma)$, and then (42) follows by the density of $\text{tr}_{\Gamma}^*(\mathbb{L}_2(\Gamma))$ in $H_{\Gamma}^{-1}$, which holds under Assumption 3.12 by Corollary 3.13. $\square$

Our method for solving the IE (20) uses $V_N$ as the approximation space in a Galerkin method, based on (21), with $a$ defined by (17). Given $g \in (\tilde{H}^1(\Gamma^c))^\perp$ we seek $\phi_N \in V_N$ such that

$$ a(\phi_N, \psi_N) = \langle g, \psi_N \rangle_{H^1(\mathbb{R}^n) \times H^{-1}(\mathbb{R}^n)}, \quad \forall \psi_N \in V_N. \quad (43) $$

Let $\{f^i\}_{i=1}^N$ be a basis for $\mathbb{V}_N$, and let $\{e^i = \text{tr}_{\Gamma}^* f^i\}_{i=1}^N$ be the corresponding basis for $V_N$. Then, writing $\phi_N = \sum_{j=1}^N c_j e^j$, (43) implies that $\vec{c} = (c_1, \dots, c_N)^T \in \mathbb{C}^N$ satisfies the system

$$ \underline{A}\vec{c} = \vec{b}, \quad (44) $$

where, by (34), (27), and (38), the matrix $\underline{A} \in \mathbb{C}^{N \times N}$ has $(i, j)$-entry given by

$$ \begin{aligned} A_{ij} &= a(e^j, e^i) = \langle \underline{A} f^j, f^i \rangle_{\mathbb{H}^{t_d}(\Gamma) \times \mathbb{H}^{-t_d}(\Gamma)} \\ &= \int_{\Gamma} \int_{\Gamma} \Phi(x, y) f^j(y) \overline{f^i(x)} \, d\mathcal{H}^d(y) \, d\mathcal{H}^d(x), \end{aligned} \quad (45) $$
and, by (30), the vector $\vec{b} \in \mathbb{C}^N$ has ith entry given by

$$b_i = \langle g, e^i \rangle_{H^1(\mathbb{R}^n) \times H^{-1}(\mathbb{R}^n)} = \int_{\Gamma} \mathrm{tr}_{\Gamma} g(x) \overline{f^i(x)} \mathrm{d}\mathcal{H}^d(x), \quad (46)$$

with $\mathrm{tr}_{\Gamma} g(x) = -u^i(x)$, $x \in \Gamma$, for the scattering problem with $g$ given by (11).

**Remark 4.2** (Connection to known cases: III). Building on Remark 3.17, if $\Gamma$ is the boundary of a bounded Lipschitz open set, or a multi-screen, then the above Galerkin method is simply a classical piecewise-constant boundary element method for the single-layer equation $S\phi = g$. If $\Gamma$ is a planar screen in the sense of Lemma 3.11 then the method is identical to that proposed in [5, §5] and the linear system (44) is identical to [5, Eqn. (55)].

Once we have computed $\phi_N$ by solving (44) we will compute approximations to $u(x)$ and $u^\infty(x)$, given by (19)/(15) and (23), respectively. Each expression takes the form $J(\phi)$, where

$$J(\psi) := \langle \varphi, \bar{\psi} \rangle_{H^1(\mathbb{R}^n) \times H^{-1}(\mathbb{R}^n)}, \quad \psi \in H_{\Gamma}^{-1}, \quad (47)$$

for some $\varphi \in (\tilde{H}^1(\Omega))^\perp$. Explicitly,

$$\varphi = P(\sigma v), \quad (48)$$

where $\sigma$ is any element of $C_{0,\Gamma}^\infty$ (with $x$ not in the support of $\sigma$ in the case that $J(\phi) = u(x)$) and $v = \Phi(x, \cdot)$ in the case that $J(\phi) = u(x)$, $v = \Phi^\infty(\hat{x}, \cdot)$ in the case that $J(\phi) = u^\infty(\hat{x})$; note that each $v$ is $C^\infty$ in a neighbourhood of $\Gamma$. In each case we approximate $J(\phi)$ by $J(\phi_N)$ which, recalling (30), is given explicitly by (cf. [5, Eqn. (62)])

$$J(\phi_N) = \langle \varphi, \bar{\phi_N} \rangle_{H^1(\mathbb{R}^n) \times H^{-1}(\mathbb{R}^n)} = \bar{c}^T \bar{\varphi}, \quad (49)$$

where $\bar{\varphi}$ has jth entry given by

$$\bar{\varphi}_j = \int_{\Gamma} \mathrm{tr}_{\Gamma} \varphi(x) f^j(x) \mathrm{d}\mathcal{H}^d(x), \quad (50)$$

and $\mathrm{tr}_{\Gamma} \varphi(x) = v(x)$, $x \in \Gamma$, for $\varphi$ given by (48). The following is a basic convergence result.

**Theorem 4.3.** Let $\Gamma$ be a compact $d$-set for some $n - 2 < d \le n$, and suppose that Assumptions 3.1 and 3.12 hold. Suppose also that $h \to 0$ as $N \to \infty$. Then for sufficiently large $N \in \mathbb{N}$ the variational problem (43) has a unique solution $\phi_N \in V_N$ that is quasi-optimal in the sense that, for some constant $C > 0$ independent of $\phi$ and $N$,

$$\|\phi - \phi_N\|_{H^{-1}(\mathbb{R}^n)} \le C \inf_{\psi_N \in V_N} \|\phi - \psi_N\|_{H^{-1}(\mathbb{R}^n)}, \quad (51)$$

where $\phi \in H_{\Gamma}^{-1}$ denotes the solution of (20). Furthermore, $\|\phi - \phi_N\|_{H^{-1}(\mathbb{R}^n)} \to 0$ as $N \to \infty$, and, where $J(\cdot)$ is given by (47) for some $\varphi \in (\tilde{H}^1(\Omega))^\perp$, $J(\phi_N) \to J(\phi)$ as $N \to \infty$.
Proof. The sesquilinear form $a(\cdot, \cdot)$ is compactly perturbed coercive (Lemma 3.3), and invertible if Assumption 3.1 holds (Theorem 3.4), so the quasi-optimality (51) holds for all sufficiently large $N$ by (42) and standard Galerkin method theory [33, §4.2.3]. The remaining results follow by Lemma 4.1 and the continuity of the linear functional $J(\cdot)$. $\square$

## 4.1 Galerkin error estimates

If the exact solution $\phi$ possesses sufficient regularity and the spaces $V_N$ have appropriate approximability properties, then Theorem 4.3 can be used to derive Galerkin error estimates, and superconvergence estimates for functionals. We record this fact in the following theorem.

**Theorem 4.4.** Let the assumptions of Theorem 4.3 hold. Suppose additionally that $\phi \in H_{\Gamma}^s$ for some $-1 < s < -(n-d)/2$, and that

$$ \inf_{\psi_h \in V_N} \|\psi - \psi_h\|_{H_{\Gamma}^{-1}} \le ch^{1+s} \|\psi\|_{H_{\Gamma}^s}, \quad 0 < h \le \text{diam}(\Gamma), \quad \psi \in H_{\Gamma}^s. \quad (52) $$

Then, for some constant $c > 0$ independent of $h$ and $\phi$,

$$ \|\phi - \phi_N\|_{H_{\Gamma}^{-1}} \le ch^{1+s} \|\phi\|_{H_{\Gamma}^s}, \quad (53) $$

for all sufficiently large $N$. Furthermore, let $J(\cdot)$ be given by (47) for some $\varphi \in (\tilde{H}^1(\Omega))^{\perp}$, and denote by $\zeta \in H_{\Gamma}^{-1}$ the solution $\phi$ of (21), in the case that $g$ is replaced by $\varphi$. Suppose that $\zeta \in H_{\Gamma}^s$. Then

$$ |J(\phi) - J(\phi_N)| \le ch^{2(1+s)} \|\phi\|_{H_{\Gamma}^s} \|\zeta\|_{H_{\Gamma}^s}, \quad (54) $$

for some constant $c > 0$ independent of $h$, $\phi$, and $\zeta$, for all sufficiently large $N$.

**Proof.** The bound (53) follows from (52) and (51). The bound (54) follows from (53) by a standard Aubin-Nitsche trick argument, as in the proof of [5, Thm 5.6]. $\square$

In the case where $\Gamma$ is the attractor of an OSC-IFS there is a natural way to build quasi-uniform meshes on $\Gamma$. Furthermore, under certain assumptions, we prove in Theorem 4.5 that the conditions of Theorem 4.4 are satisfied, so that the error bounds (53) and (54) hold.

Let $\Gamma \subset \mathbb{R}^n$ be the attractor of an OSC-IFS $\{s_1, \dots, s_M\}$. Following [23], for $\ell \in \mathbb{N}$ we define the set of multi-indices $I_\ell := \{1, \dots, M\}^\ell = \{\mathbf{m} = (m_1, m_2, \dots, m_\ell)\}$: $1 \le m_j \le M, j = 1, 2, \dots, \ell$}, and for $E \subset \mathbb{R}^n$ and $\mathbf{m} \in I_\ell$ we define $E_\mathbf{m} = s_{m_1} \circ s_{m_2} \circ \dots \circ s_{m_\ell}(E)$. We also set $I_0 := \{0\}$ and adopt the convention that $E_0 := E$. This notation extends that of (4) where the sets $\Gamma_1, \dots, \Gamma_M$ were introduced, corresponding to $E = \Gamma$ and $\ell = 1$.
<figure><img src="image_5.png" /><figcaption>(a) $\ell = 1$ decomposition</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) $\ell = 2$ decomposition</figcaption></figure>

Figure 4: Level 1 (a) and level 2 (b) decompositions of the Koch curve. To make the labelling more compact, in (a) the labels “1”, . . . , “4” indicate the subsets Γ$_{1}$, . . . , Γ$_{4}$, and in (b) the labels “ij” and “ijk” indicate Γ$_{(}$$_{i,j}$$_{)}$ and Γ$_{(}$$_{i,j,k}$$_{)}$ . In (b) the insert shows the level 3 decomposition of Γ$_{(1}$$_{,}$$_{1)}$ .

here. We illustrate this for the Koch curve (Example 2.3, Figure 1(e)) in
Figure 4. Illustrations for other examples are given in [17, Figs 1-6].

Let $0 < h \le \operatorname{diam}(\Gamma)$. Define the index set $L_h$ by $L_h := \{0\}$ for $h = \operatorname{diam}(\Gamma)$, and by

$$
L_h := \left\{ \boldsymbol{m} \in \bigcup_{\ell=1}^{\infty} I_{\ell} : \operatorname{diam}(\Gamma_{\boldsymbol{m}}) \le h \text{ and } \operatorname{diam}(\Gamma_{\boldsymbol{m}_{-}}) > h \right\}, \quad (55)
$$

for $h < \operatorname{diam}(\Gamma)$, where, for $\boldsymbol{m} = (m_1, \dots, m_\ell)$, $\boldsymbol{m}_- := (m_1, \dots, m_{\ell-1})$ if
$\ell \in \mathbb{N}$ with $\ell \ge 2$, and $\boldsymbol{m}_- := 0$ if $\ell = 1$. Then $\{T_j\}_{j=1}^N := \{\Gamma_m\}_{m \in L_h}$
defines a quasi-uniform mesh of $\Gamma$. We define the spaces of piecewise-
constant functions

$$
\mathbb{Y}_h := \operatorname{span}(\{\chi_m\}_{m \in L_h}) \quad \text{and} \quad Y_h := \operatorname{tr}^*_{\Gamma}(\mathbb{Y}_h) \subset H_{\Gamma}^{-1}, \quad (56)
$$

where {$\chi_m$}$_{m \in L_h}$ is the canonical $\mathbb{L}_2(\Gamma)$-orthonormal basis for $\mathbb{Y}_h$ given by

$$
\chi_m(x) := \begin{cases} \frac{1}{H^d(\Gamma_m)^{1/2}}, & x \in \Gamma_m, \\ 0, & \text{otherwise,} \end{cases} \tag{57}
$$

The following theorem is then a consequence of results in [5] and [4].

**Theorem 4.5.** Let $\Gamma$ be an OSC-IFS attractor with dimension $d := \dim_{\mathbb{H}}(\Gamma)$, such that either (a) $\Gamma$ is disjoint, (b) $d = n$, or (c) $d = n-1$ and $\Gamma \subset \Gamma_\infty := \mathbb{R}^{n-1} \times \{0\}$. Let Assumption 3.1 hold, and set $V_N = Y_h$. Then (52) holds for all $-1 < s < -(n-d)/2$.

Furthermore, for the scattering problem with g defined by (11) there ex-
ists $-1 < s < -(n-d)/2$ such that (53) holds and (54) holds for both
$J(\phi) = u(x)$ and $J(\phi) = u^\infty(\hat{x})$.
*Proof.* For case (a), where $\Gamma$ is a disjoint IFS attractor, it was proved in [5, Prop. 5.2] that, for every $0 < t < 1$ and every $0 < t' < t$, there exists a constant $c > 0$ such that

$$ \inf_{\Psi_h \in \mathbb{Y}_h} \|f - \Psi_h\|_{\mathbb{H}^{-t}(\Gamma)} \le c h^{t-t'} \|f\|_{\mathbb{H}^{-t'}(\Gamma)}, \quad 0 < h \le \mathrm{diam}(\Gamma), \quad f \in \mathbb{H}^{-t'}(\Gamma). \tag{58} $$

In [5, Prop. 5.2] this result was actually only stated for $n-1 < d < n$, but the argument of [5, Prop. 5.2] holds in fact for all $0 < d < n$ (so, in particular for $n-2 < d < n$), because the results from [23] on which it is based hold for all $0 < d < n$. This latter statement requires some explanation. A key step in the argument of [5, Prop. 5.2] was the use of results from [23] to prove that $\|\cdot\|_{\mathbb{H}^t(\Gamma)}$ is equivalent to a norm defined in terms of coefficient decay in a wavelet expansion; see [5, Thm. 3.1 and Cor. 3.3(iii)]. The relevant results in [23] (Theorems 1 and 2) are stated under the additional assumption that $\Gamma$ is not contained in an $(n-1)$-dimensional hyperplane. However, this additional assumption is made in [23] solely to ensure that Markov's inequality [23, Eqn (4.1)] is satisfied for whatever class of polynomials is being used in the wavelet expansion. As a result, this additional assumption is superfluous for us because we consider only piecewise-constant approximations and $0 < t < 1$, while [23] considers also higher order polynomials and larger $t$, and for constant functions Markov's inequality [23, Eqn (4.1)] is trivially satisfied.

For case (b), the result (58) was proved for $t=1$ in [4], using a quite different argument based on Poincaré inequalities. The fact that it also holds in case (c), again for $t=1$, follows from the result for case (b), applied in the setting of $\mathbb{R}^{n-1}$, and Lemma 3.11.

The above establishes (58) for the particular case $t = t_d$. The bound (52) then follows by Theorem 3.9 and Corollary 3.13, noting that in case (b) Assumption 3.12 holds by Remark 3.14(iii).

The final statement then follows from Theorem 4.4, since, for the scattering problem, and the choices of $J$ under consideration, the solutions $\phi$ and $\zeta$ possess some extra regularity by Proposition 3.19 (see the argument in Remark 3.20 for $\phi$, and argue similarly for $\zeta$). $\square$

**Remark 4.6 (Convergence rates).** Suppose that, in addition to the assumptions of Theorem 4.5, Hypothesis 3.21 holds. Then, arguing as in Remark 3.20, assuming the datum $g$ is sufficiently smooth, we will have the maximum possible regularity for $\phi$, i.e., $\phi \in H_\Gamma^s$ for every $-1 < s < -(n-d')/2$. Then, assuming that the bounds in Theorem 4.4 are sharp, in numerical experiments we should expect to see errors in the computation of $\phi$ and of linear functionals of $\phi$ roughly proportional to $h^{1+(d'-n)/2}$ and $h^{2+d'-n}$, respectively. (For the latter, assume also that $\varphi$ in Theorem 4.3 is sufficiently smooth so that $\zeta \in H_\Gamma^s$ for every $-1 < s < -(n-d')/2$.)

If, additionally, $\Gamma$ is homogeneous, with $\rho_m = \rho$ for $m=1, \dots, M$, for some $0 < \rho < 1$, in which case $d = \dim_{\mathbb{H}}(\Gamma) = \log(M)/\log(1/\rho)$, then
the meshes defined by (55) are uniform, and taking $h = \rho^\ell$ diam $\Gamma$, for some $\ell \in \mathbb{N}$, gives $L_h = I_\ell = \{1, \dots, M\}^\ell$ and $V_N = Y_h = \text{tr}_\Gamma^*(\text{span}(\{\chi_m\}_{m \in I_\ell}))$. In this case, since $h$ is proportional to $\rho^\ell$ and $\rho = M^{-1/d}$, we should see errors in $\phi$ and in linear functionals of $\phi$ proportional to $(M^{d'/d})^{-\ell/2}$ and $(M^{d'/d})^{-\ell}$, respectively, in the case $n=2$, and proportional to $(M^{d'/d}\rho)^{-\ell/2}$ and $(M^{d'/d}\rho)^{-\ell}$, respectively, in the case $n=3$.

## 4.2 Numerical quadrature

To implement our method we need suitable numerical quadrature rules to evaluate the integrals (45), (46) and (50). For this we generalise the approach taken for the screen case in [5]. Here we give only the main ideas, and refer the reader to Appendix A, [5, §5.4], and [18, 17] for details.

Suppose that $\Gamma$ is an OSC-IFS attractor, and that, as in §4.1, we are using the approximation space $V_N = Y_h$ given by (56). Suppose that $g$ is given by (11) and $\varphi$ by (48), with $u^i$ and $v$ both $C^\infty$ in a neighbourhood of Hull($\Gamma$), the convex hull of $\Gamma$. Suppose that we adopt the canonical $L_2(\Gamma)$-normalised basis (57), so that $f^j = \chi_{m(j)}$, $j = 1, \dots, N$, where $N := |L_h|$, with $L_h$ given by (55), and $(m(1), \dots, m(N))$ is some ordering of the elements of $L_h$. Then, where $\mu_m := \mathcal{H}^d(\Gamma_m)$ for $m \in L_h$, the integrals to be evaluated are, for $i, j \in \{1, \dots, N\}$,

$$A_{ij} = \mu_{m(i)}^{-1/2} \mu_{m(j)}^{-1/2} \int_{\Gamma_{m(i)}} \int_{\Gamma_{m(j)}} \Phi(x, y) \, d\mathcal{H}^d(y) \, d\mathcal{H}^d(x), \quad (59)$$

$$b_i = -\mu_{m(i)}^{-1/2} \int_{\Gamma_{m(i)}} u^i(x) \, d\mathcal{H}^d(x), \qquad \vec{\varphi}_i = \mu_{m(i)}^{-1/2} \int_{\Gamma_{m(i)}} v(x) \, d\mathcal{H}^d(x). \qquad (60)$$

Since $u^i$ and $v$ are smooth in a neighbourhood of $\Gamma$, (60) can be evaluated using the composite barycentre rule of [18, Defn 3.1], cf. [5, (97)-(99)]. This involves decomposing the mesh element $\Gamma_{m(i)}$ into smaller self-similar sub-elements whose vector indices are taken from the index set $L_{h_Q}$, for some maximum quadrature element diameter $h_Q \le h$, and applying a one-point quadrature rule on each sub-element. Similarly, provided that $\Gamma_{m(i)}$ and $\Gamma_{m(j)}$ are disjoint, (59) can be evaluated using a tensor product version of this composite barycentre rule (defined in [18, Defn 3.5]), cf. [5, (92)].

When $\Gamma_{m(i)}$ and $\Gamma_{m(j)}$ are not disjoint, the integral in (59) is singular. Singularity subtraction reduces the problem to the evaluation of

$$\int_{\Gamma_{m(i)}} \int_{\Gamma_{m(j)}} \Phi_{\text{sing}}(x, y) \, d\mathcal{H}^d(y) \, d\mathcal{H}^d(x), \quad (61)$$

where $\Phi_{\text{sing}}(x, y) = -\log(|x - y|)/(2\pi)$ if $n=2$, and $\Phi_{\text{sing}}(x, y) = 1/(4\pi|x - y|)$ if $n=3$. The integral of $\Phi - \Phi_{\text{sing}}$ is regular and can be evaluated using
the tensor product composite barycentre rule, cf. [5, (94)]. The treatment of (61) depends on the nature of $\Gamma$.

If $\Gamma$ is disjoint (e.g. the Cantor set, Figure 1(e)) then (61) is singular if and only if $i = j$, in which case (61) can be evaluated using the quadrature rules of [18, §4.3], cf. [5, (95)-(96)]. These rules exploit the self-similarity of $\Gamma$ and the homogeneity of $\Phi_{\text{sing}}$ to write the singular integral (61) in terms of regular integrals, which can be evaluated by the composite barycentre rule.

If $\Gamma$ is non-disjoint (e.g. the Koch curve, or the Koch snowflake, Fig. 1(f), (g)) then the situation is more complicated, because, in addition to the self-interaction case $i = j$, (61) can also be singular for $i \neq j$, if $\Gamma_{m(i)}$ and $\Gamma_{m(j)}$ intersect at a point or at a higher-dimensional set. For certain non-disjoint attractors, it holds that: (i) all singular instances of (61) that arise in our discretization can be written in terms of one of a finite collection of “fundamental” singular integrals, which capture the different singular interactions that can occur between mesh elements; and (ii) these fundamental singular integrals together satisfy a small linear system of equations that can be solved in closed form in terms of regular integrals, which can be evaluated using the composite barycentre rule. A general algorithm for identifying the fundamental singular integrals and deriving the associated linear system was presented in [17, Algorithm 1], along with explicit formulas for the Sierpinski triangle, the Vicsek fractal, the Sierpinski carpet, and the Koch snowflake. These formulas were applied in the context of screen scattering problems in [17, §7.3]. In Appendix A we briefly explain the methodology of [17], and derive explicit formulas for the case of the Koch curve, which was not considered in [17].

The accuracy of the quadrature approximations described above for the evaluation of (59) and (60) can be controlled by a single parameter $h_Q \le h$, which represents the maximum diameter of the sub-elements used in the composite barycentre rule. Using the results of [18] one can prove quadrature error estimates. The following theorem is a generalisation of [5, Thm 5.14]. While [5, Thm 5.14] was stated for the special case where $\Gamma \subset \mathbb{R}^{n-1} \times \{0\}$ is a planar screen, it extends trivially to our more general context, with minor notational adjustments, because the planarity of $\Gamma$ was not used in its proof. We recall that $\text{Hull}(E)$ denotes the convex hull of $E \subset \mathbb{R}^n$, and we denote by $\|\cdot\|_2$ both the Euclidean norm on $\mathbb{C}^N$ and the induced matrix norm on $\mathbb{C}^{N \times N}$.

**Theorem 4.7.** Let $\Gamma$ be an OSC-IFS attractor. Let $A^Q$, $\tilde{b}^Q$ and $\tilde{\varphi}^Q$ denote the approximations of (59) and (60) obtained via the quadrature described above, using a maximum sub-element diameter $0 < h_Q \le h$ in the composite barycentre rule.

(i) Let $u^i$ satisfy the Helmholtz equation in some open neighbourhood of
Hull(Γ). Then

$$ \|\vec{b} - \vec{b}^Q\|_2 \le h_Q^2 |u^i|_{2, \text{Hull}(\Gamma)} \mathcal{H}^d(\Gamma)^{1/2}, \quad (62) $$

where $|u^i|_{2, \text{Hull}(\Gamma)} := \max_{x \in \text{Hull}(\Gamma)} \max_{\substack{\alpha \in \mathbb{N}_0^N \\ |\alpha|=2}} |D^\alpha u^i(x)|$.

(ii) Let $v$ be $C^\infty$ in a neighbourhood of Hull(Γ). Given $\psi_N \in V_N = Y_h$, let $J^Q(\psi_N)$ be defined by (49) with $\vec{\phi}$ replaced by $\vec{\phi}^Q$, and let $\vec{\psi}$ denote the coefficient vector of $\psi_N$ in the basis $\{e^i = \text{tr}_\Gamma^* f^i\}_{i=1}^N$. Then, there exists $C > 0$, independent of $h$, $h_Q$, $v$ and $\psi_N$, such that

$$ |J(\psi_N) - J^Q(\psi_N)| \le h_Q^2 |v|_{2, \text{Hull}(\Gamma)} \|\vec{\psi}\|_2 \mathcal{H}^d(\Gamma)^{1/2}. \quad (63) $$

(iii) Suppose that $\Gamma$ is hull-disjoint, meaning that $\text{Hull}(\Gamma_m) \cap \text{Hull}(\Gamma_{m'}) = \emptyset$ for every $m \ne m' \in \{1, \dots, M\}$. Then there exists $C > 0$, independent of $h$ and $h_Q$, such that

$$ \|\underline{A} - \underline{A}^Q\|_2 \le Ch_Q h^{-(n-1)} \mathcal{H}^d(\Gamma). \quad (64) $$

If, further, $\Gamma$ is homogeneous, then

$$ \|\underline{A} - \underline{A}^Q\|_2 \le Ch_Q^2 h^{-n} \mathcal{H}^d(\Gamma). \quad (65) $$

While Theorem 4.7(iii) is stated only for hull-disjoint attractors (because that was the setting considered in [5]), we expect it should be possible to prove similar results for non-disjoint OSC-IFS attractors, by combining the results of [18, §4.3] with those of [17, §6]. Indeed, numerical experiments (not reported here) suggest that $\|\underline{A} - \underline{A}^Q\|_2 = O(h_Q^2)$ as $h_Q \to 0$ for all the examples considered in §5, including the Koch snowflake, which is both non-disjoint and non-homogeneous. However, we leave the proof of this for future work.

In principle, the quadrature error estimates of Theorem 4.7 could be combined with the semi-discrete convergence estimates of Theorem 4.4 to obtain a fully discrete analysis for our IE method (under appropriate assumptions, such as disjointness), with conditions on how small $h_Q$ should be in order to maintain the convergence rates in Theorem 4.4. For brevity we do not embark on such an analysis here, but refer the interested reader to [5, §5.4] where the analogous analysis was carried out for screen problems. In practice, our numerical results in §5 suggest that in many cases it may be sufficient to decrease $h_Q$ in proportion to $h$ in order to achieve the predicted rates.
# 5 Numerical results

In this section we present numerical results obtained using our Galerkin IEM for scattering by various fractals $\Gamma$, each an OSC-IFS attractor$^6$. For each example we assume plane wave incidence, i.e. the datum $g$ is as in (11) with $u^i(x) = e^{ik\vartheta \cdot x}$ and $|\vartheta| = 1$, and we compute the Galerkin IEM solution by solving (44), using the piecewise-constant quasi-uniform-mesh approximation space $V_N = Y_h$, with $Y_h$ defined as in (56), so that each element is a scaled copy of $\Gamma$, and with the basis functions as defined above (59). We approximate the scattered field $u$ and/or the far field $u^\infty$, which each (as discussed above (47)) take the form of the linear functional (47) with $\varphi$ given by (48), by the discretisation (49). These calculations require evaluation of the integrals (45), (46), and (50). We approximate these by the methods detailed in §4.2, using a maximum quadrature element size $h_Q = C_Q h$. We choose $C_Q := \max\{\rho_m^2 : m = 1, \dots, M\}$, where $\rho_m$, $m = 1, \dots, M$, are the contraction factors of the IFS, except for the higher wavenumber simulations for $k \ge 20$ in Figures 5, 7 and 8, where we use $C_Q := \max\{\rho_m^4 : m = 1, \dots, M\}$. To validate the accuracy of our quadrature, a number of our experiments were repeated using smaller values for $h_Q$, and the difference in the results was found to be negligible.

When we plot errors we use as our “exact” solution a more accurate Galerkin-IEM solution. Most of our experiments are for homogeneous attractors, in which case our mesh is uniform with $N = M^\ell$, for some $\ell \in \mathbb{N}$, we denote the corresponding approximate scattered, total, and far fields by $u_\ell$, $u_\ell^t := u^i + u_\ell$, and $u_\ell^\infty$, respectively, and the “exact” solution is the solution for $\ell = \ell_{\text{ref}}$, for some $\ell_{\text{ref}}$ that we note for each example. Where we plot $L^\infty$ relative error estimates these are

$$ \frac{\|u_{\ell_{\text{ref}}} - u_\ell\|_{L^\infty}}{\|u_{\ell_{\text{ref}}}\|_{L^\infty}} \quad \text{and} \quad \frac{\|u_{\ell_{\text{ref}}}^\infty - u_\ell^\infty\|_{L^\infty}}{\|u_{\ell_{\text{ref}}}^\infty\|_{L^\infty}}, \qquad (66) $$

where the $L^\infty$ norms are discrete norms taken over a set of points detailed for each example.

## 5.1 Examples in 2D space

### 5.1.1 Cantor dust and Koch curve

Plots of Re($u_\ell^t$) are shown in Figure 5 for two fractal scatterers. The first, see (a)-(c), is the middle-third Cantor dust, $\Gamma = C \times C$, where C is the Cantor set defined in Example 2.2, and the other, panel (d), is the Koch curve of Example 2.3. For both scatterers the IFS is homogeneous, with $M = 4$, $\rho = 1/3$, and hence $d = \dim_H(\Gamma) = \log 4 / \log 3 \approx 1.26$; see, e.g., [5,

<footnote>$^6$Our method is implemented in Julia and is available to download at github.com/AndrewGibbs/IFSIntegrals</footnote>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: Scattering by a middle third Cantor dust, (a)-(c), and a Koch curve, (d). See §5.1(5.1.1).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6: Plots of the discrete relative-error estimates (66) for a range of 2D examples; see §5.1(5.1.2).</figcaption></figure>

Eqn. (125)] for the Cantor dust IFS. In all plots the incident plane wave has direction $\vartheta = (1, 1)/\sqrt{2}$, and we take $k = 20$ and $\ell = 4$ in (a) and (d), $k = 60$ and $\ell = 5$ in (b)-(c), so that $h = 1/3^4 \approx 0.0123$ in (d), $h = \sqrt{2}/3^4 \approx 0.0175$ in (a), and $h = \sqrt{2}/3^5 \approx 0.00582$ in (b)-(c). The key difference between (a) and (b) is the tripling of $k$, so that the wavelength $\lambda = 2\pi/k$ reduces from $\lambda \approx 0.314$ in (a) to $\lambda \approx 0.105$ in (b). The wave field in (a) does not appear to resolve details beyond level 2, i.e. it appears that $u^t \approx 0$ in the convex hull of each of the sixteen $\Gamma_m$ with $m \in I_2$ (in the notation of §4.1). This is unsurprising as $u^t = 0$ on each level 2 component, $\Gamma_m$, with $m \in I_2$, and each is comprised of four level 3 components on which $u^t = 0$ and whose separation is only $1/3^3 \approx 0.118\lambda$, i.e., is a small fraction of $\lambda$. In (b), where $\lambda$ is reduced by a factor 3, the wave field appears to resolve detail down to level 3, i.e. to resolve details of $1/3$ the size. To see this more clearly the region inside the dotted boundary is blown up by a factor 3 in (c). After this scaling in fact, thanks to the incidence direction we have chosen, the part of the plot (c) in $[-1/3, 1/3]^2$ is very similar to the field plotted in (a) in $[-1, 1]^2$.

### 5.1.2 Convergence plots

In Figure 6 we show the discrete $L^\infty$ relative errors (66) for a range of 2D examples, namely the Koch curve (Ex. 2.3), the Cantor set $C \times \{0\}$ (Ex. 2.2), and the Cantor dust $C \times C$ with two different values of $\rho$, with plane wave incidence direction $\vartheta = (1, -1)/\sqrt{2}$ and wavenumber $k = 5$. To compute the scattered-field relative error given by (66) we sample at 50 points equispaced along each edge of the square $(-1, 2) \times (-1.5, 1.5)$ (200 points in total), and for the far-field we sample at 50 equispaced points on the circle $S^1$. We use, for each scatterer, $\ell_{\text{ref}} = \ell_{\text{max}} + 2$, where $\ell_{\text{max}}$ is the largest $\ell$ for which
results are shown.

Also plotted in Figure 6 are graphs of $cM^{-\ell}$, with $c > 0$ chosen to fit each error curve. In the cases where Theorem 4.5 applies (all except the Koch curve), then if Hypothesis 3.21 holds we expect, by Remark 4.6, errors to be roughly proportional to $M^{-\ell}$. In the examples with $d \approx 0.631$ the relative errors do seem to be proportional to $M^{-\ell}$ for larger $\ell$, supporting Hypothesis 3.21 in these cases. The errors in the two examples with $d \approx 1.262$ seem to decrease at the same rate, suggesting the same solution regularity in both cases, and that the error estimate of Theorem 4.5 may hold also for the Koch curve even though $\Gamma$ is not disjoint in that case. But the convergence is slower than $M^{-\ell}$, suggesting that Hypothesis 3.21 does not hold in these cases.

### 5.1.3 The Koch snowflake

In Figure 7 we show approximations to $\text{Re}(u^t)$ for scattering by a Koch snowflake for the same incident plane wave as Figure 5, computed in two different ways, illustrating Remark 3.8. In Figure 7(a) we solve the IE by our Galerkin IEM with $N = 4039$ on the solid Koch snowflake $\Gamma$, shown in Fig. 1(g), which is the attractor of a non-homogeneous IFS with $M = 7$ as noted in Ex. 2.4. We refer to this as the *volume approach*. In Figure 7(b) we solve the IE by our Galerkin IEM on $\partial\Gamma$, the boundary of the snowflake. We refer to this as the *boundary approach*. In contrast to all our other examples, $\partial\Gamma$ is not an IFS attractor, but it is the union of three IFS attractors (rotated copies of the Koch curve of Ex. 2.3, each the attractor of an IFS with $M = 4$), and so $\partial\Gamma$ is a $d'$-set, with $d' := \dim_H(\partial\Gamma) = \log(4)/\log(3) \approx 1.262$. In Figure 7(b) we use $M^4 = 256$ degrees of freedom on each Koch curve comprising $\partial\Gamma$, so that $N = 768$.

In the boundary approach, to assemble the Galerkin matrix $\underline{A}$, we view it as a $3 \times 3$ block matrix, each block corresponding to interactions between two of the three Koch curves. The diagonal blocks correspond to self-interactions for a single Koch curve, and these matrix elements are approximated by quadrature as described at the beginning of the section. The off-diagonal blocks are assembled using the composite barycentre rule using the same value of $h_Q$ as for the diagonal blocks.

Proposition 3.6 and Remark 3.2 tell us that the IE solution in the volume approach is supported on $\partial\Gamma$, and that the IE solutions and scattered fields for the two approaches coincide as long as $k^2$ is not a Dirichlet eigenvalue of $-\Delta$ in $\Omega_- := \Gamma^\circ$, the interior of the snowflake. It appears that $k = 20$ is not one of these resonant wavenumbers as the fields in Figure 7(a) and (b) coincide and the field in $\Omega_-$ is zero in (b), in agreement with the boundary condition for the volume approach that $u^t \in \tilde{H}^1(\Omega)$, where $\Omega = \Gamma^c$.

In Figure 7(c)-(e) we plot the modulus of the piecewise-constant Galerkin IEM solution $\phi_N \in V_N = Y_h \subset H_\Gamma^{-1}$ corresponding to Figure 7(a), for
<figure><img src="image_10.png" /><figcaption>(a) Re($u^t$) for volume approach</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Re($u^t$) for boundary approach</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) $|\phi_N|$ for volume approach, $h = 0.22$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d) $|\phi_N|$ for volume approach, $h = 0.074$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(e) $|\phi_N|$ for volume approach, $h = 0.025$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(f) Far-field relative $L^\infty$ errors</figcaption></figure>

Figure 7: Scattering by the Koch snowflake; see §5.1(5.1.3). (a) and (b) show Re($u^t$) computed with both approaches, for $k = 20$. For the same $k$, (c)-(e) show $|\phi_N|$, where $\phi_N$ is the Galerkin IEM solution for the volume approach, for three different $h$ values. (f) plots the far-field relative $L^\infty$ error for both approaches against the total number of degrees of freedom, for $k = 5$.
$h \approx 0.22$, $h \approx 0.074$, and $h \approx 0.025$. Since $\phi_N$ is constant on each element, the meshes used for each $h$ are discernible in Figure 7(c)-(e) (each element is a scaled copy of the original snowflake $\Gamma$). Each solution $\phi_N$ is highly peaked near $\partial\Gamma$, especially where $\partial\Gamma$ is illuminated by the incident wave, and is much smaller away from $\partial\Gamma$; these effects are increasingly marked as $h$ is reduced. This is unsurprising as, by Theorem 4.3, $\phi_N \to \phi$ in $H_\Gamma^{-1}$ as $h \to 0$, and $\phi$ is supported in $\partial\Gamma$.

In Figure 7(f) we explore convergence of the far-field approximations $u_\ell^\infty$ computed by the volume and boundary approaches, showing computations for $h = \text{diam}(\Gamma)/3^{\ell/2}$ for $\ell = 0, 1, \dots, 6$ for the volume approach, and $h = 3\text{diam}(\Gamma)/3^\ell$ for $\ell = 0, 1, \dots, 4$ for the boundary approach. For each approach we use as our “exact” solution the boundary approach solution with $\ell = \ell_{\text{ref}} = 7$. Figure 7(f) shows the relative errors (66) in $u_\ell^\infty$ for both methods, for plane wave incidence direction $\vartheta = (1, -1)/\sqrt{2}$ and $k = 5$, with the $L^\infty$ norms computed using the same discrete set of points as in Figure 6. In the volume approach every second increment in $\ell$ has a smaller reduction in error. At these increments, the elements adjacent to $\partial\Gamma$, which is the support of the solution, are not being subdivided, as a consequence of the definition of the approximation space $Y_h$. The convergence rate results of Theorem 4.5 apply to the volume approach but not to the boundary approach as $\partial\Gamma$ is not an IFS attractor. But, assuming these estimates apply in both cases, and if Hypothesis 3.21 holds, so that the solution $\phi \in H_{\partial\Gamma}^{-1}$ has its maximum possible regularity, then as in Remark 4.6 we anticipate errors decreasing roughly in proportion to $h^{d'}$ in both cases, i.e. proportional to $N^{-1}$ and $N^{-d'/2} \approx N^{-0.631}$ in the respective boundary and volume cases. Both approaches appear to be converging somewhat more slowly than these conjectured theoretical rates, but, of the two, the boundary approach is certainly converging more rapidly. It is plausible that the boundary approach, in which only $\partial\Gamma$ is discretised, should be more efficient, given that the solution is supported on $\partial\Gamma$. But the IE on $\partial\Gamma$ is not well-posed for all $k > 0$, in contrast to the IE on $\Gamma$, and there must be scope to improve the efficiency of the volume approach by using graded versions of our meshes, concentrating elements near $\partial\Gamma$ (cf. [25, 6]).

## 5.2 Examples in 3D space

In Figure 8 we show the real parts of the scattered fields created by the two Sierpinski tetradehra of Figure 2, which are attractors of the homogeneous IFS of Example 2.5 with $M = 4$ and $d = \log 4 / \log(1/\rho)$. The plane wave incidence direction is $\vartheta = (0, 1, -1) / \sqrt{2}$, $k = 50$, and both approximations were computed with $\ell = 7$, corresponding to $N = 16384$.

In Figure 9 we show $L^\infty$ far-field errors for the same incidence direction and $k = 2$ for a range of 3D examples, namely: the Sierpinski tetrahedra of Figures 2 and 8 and Example 2.5; 3D Cantor Dusts, i.e., $C \times C \times C$ where
<figure><img src="image_16.png" /><figcaption>Figure 8: Re($u_ℓ$) for scattering by the Sierpinski tetrahedra with $k = 50$ and $ℓ = 7$. See §5.2.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 9: $L^∞$ convergence plots for a range of 3D examples. See §5.2.</figcaption></figure>
C is the Cantor set of Example 2.2, with $\rho = 1/3$ and $\rho = 0.45$; the Koch curve of Figure 4 embedded in 3D space, i.e. $K \times \{0\}$, where $K \subset \mathbb{R}^2$ is the Koch curve of Example 2.3. All these scatterers $\Gamma$ have $d = \dim_H(\Gamma) > 1$ (see Figure 9) so that $H_\Gamma^{-1}$ is non-trivial by Remark 3.5, the Galerkin IEM is applicable, and the solution $\phi \in H_\Gamma^{-1}$ to the IE (20) is non-zero (since (20) is equivalent to (35) and $\mathrm{tr}_\Gamma g = -u^i|_\Gamma$ is non-zero), so also (by (14)) the scattered field $u = A\phi$ is non-zero.

To compute the discrete $L^\infty$ relative errors (66) shown in Figure 9(a) we sample at 200 points on the sphere $S^2$, chosen so that the points form a uniform grid in spherical coordinate space $[0, \pi] \times [0, 2\pi]$, and we use, for each scatterer, $\ell_{\mathrm{ref}} = \ell_{\mathrm{max}} + 1$, where $\ell_{\mathrm{max}}$ is the largest $\ell$ for which results are shown. This choice of $\ell_{\mathrm{ref}}$, constrained by computational resources, is not large enough for $u_{\ell_{\mathrm{ref}}}^\infty$ to be a sufficiently accurate “exact” solution to see convergence rates clearly. Thus we also plot in Figure 9(b) the absolute increment errors $\|u_\ell^\infty - u_{\ell+1}^\infty\|_{L^\infty}$ for $\ell = 0, \dots, \ell_{\mathrm{ref}} - 1$. As discussed in [5, §6.2], if, for some $c > 0$ and $0 < \alpha < 1$, $\|u_\ell^\infty - u_{\ell+1}^\infty\|_{L^\infty} = c\alpha^\ell$ for all $\ell \ge \ell_0$ then, by the triangle rule, $\|u^\infty - u_\ell^\infty\|_{L^\infty} \le \frac{c}{1-\alpha}\alpha^\ell$ for $\ell \ge \ell_0$. Thus convergence rates can be deduced from Figure 9(b). By Remark 4.6, which applies to all the examples except the Koch curve and the Sierpinski tetrahedron with $d = 2$, we expect, if Hypothesis 3.21 holds, to see errors roughly proportional to $(M\rho)^{-\ell}$. This rate is observed in Figure 9(b) for sufficiently large $\ell$ for all the cases with $d < 2$, but the convergence is significantly slower than $(M\rho)^{-\ell}$ for the example with $d \approx 2.6$. These results, and the convergence results reported in §5.1, suggest that Hypothesis 3.21 does not hold in cases where $d' = \dim_H(\partial\Gamma) > n - 1$ (note $\partial\Gamma = \Gamma$ for the scatterers in Figures 6 and 9), but may hold in cases where $d' < n - 1$. They suggest moreoever that Hypothesis 3.21 and the estimates of Theorem 4.4 may hold for the Koch curve screen in 3D, even though $\Gamma$ is non-disjoint in this case.

## A Singular quadrature on fractals

In this appendix we briefly outline the methodology of [18, 17] for the derivation of representation formulas for singular integrals on fractals in terms of regular integrals. We assume throughout that $\Gamma$ is the attractor of an IFS satisfying the OSC.

The basic singular integral we consider is

$$I_{\Gamma, \Gamma} := \iint_{\Gamma} \tilde{\Phi}_t(|x - y|) d\mathcal{H}^d(x)d\mathcal{H}^d(y),$$

where

$$\tilde{\Phi}_t(r) := \begin{cases} r^{-t}, & t > 0, \\ \log r, & t = 0. \end{cases}$$
The cases $t=0$ and $t=1$ are those relevant to §4(b), since $\Phi_{\text{sing}}(x, y) = -\tilde{\Phi}_0(|x-y|)/(2\pi)$ for $n=2$, and $\Phi_{\text{sing}}(x, y) = \tilde{\Phi}_1(|x-y|)/(4\pi)$ for $n=3$, but for completeness we consider the general case. The integral $I_{\Gamma,\Gamma}$ is finite if and only if $t < d$, where $d$ is the Haudsorff dimension of $\Gamma$ (see, e.g., [18, Cor. A.2] or [5, Cor. 2.3]), so we assume henceforth that $0 \le t < d$.

## A.1 Similarity

The representation formulas and quadrature rules presented in [18, 17] are based on decomposing $I_{\Gamma,\Gamma}$ as a sum of integrals over self-similar subsets of $\Gamma$. For $m, m' \in \bigcup_{\ell=0}^{\infty} I_{\ell}$ we define

$$I_{m,m'} := \int_{\Gamma_m} \int_{\Gamma_{m'}} \tilde{\Phi}_t(|x-y|) d\mathcal{H}^d(x)d\mathcal{H}^d(y),$$

which is singular when $\Gamma_m \cap \Gamma_{m'}$ is non-empty, and regular otherwise. Key to the methodology of [18, 17] is that many of these integrals $I_{m,m'}$, for different choices of $m, m'$, can be related to each other using the self-similarity of $\Gamma$ and the homogeneity of $\tilde{\Phi}_t$, namely that, for $\rho > 0$,

$$\tilde{\Phi}_t(\rho r) = \begin{cases} \rho^{-t} \tilde{\Phi}_t(r), & t > 0, \\ \log \rho + \tilde{\Phi}_t(r), & t = 0. \end{cases} \qquad (67)$$

The following result is a consequence of [17, Props 3.2 & 3.3, Rem. 2.1]. Here $|m|$ denotes the length of the vector index $m$, with $|0|$ interpreted as 0 in the case $\Gamma_m = \Gamma_0 = \Gamma$, and we define $\vartheta_t$ by $\vartheta_t := 0$ for $t > 0$ and $\vartheta_t := \mathcal{H}^d(\Gamma)^2$ for $t = 0$. Adopting the terminology of [17], when the conditions of Proposition A.1 hold we say that the integrals $I_{m,m'}$ and $I_{n,n'}$ are similar. The condition (68) stipulates that there exist similarities mapping $\Gamma_m$ to $\Gamma_n$, and $\Gamma_{m'}$ to $\Gamma_{n'}$, defined in terms of the IFS and symmetry properties of $\Gamma$, that are compatible in a certain sense. We note that one can always take the isometries $T$ and $T'$ to be the identity map in the following.

**Proposition A.1.** Let $m, m', n, n' \in \bigcup_{\ell=0}^{\infty} I_{\ell}$. Let $T$ and $T'$ be isometries of $\mathbb{R}^n$ such that $T(\Gamma) = T'(\Gamma) = \Gamma$. Suppose there exists $\varrho > 0$ such that

$$|s_m(T(s_n^{-1}(x))) - s_{m'}(T'(s_{n'}^{-1}(y)))| = \varrho|x-y|, \quad x, y \in \mathbb{R}^n. \qquad (68)$$

Then

$$\varrho = \frac{\rho_m}{\rho_n} = \frac{\rho_{m'}}{\rho_{n'}} \qquad (69)$$

and

$$I_{m,m'} = \varrho^{2d-t} I_{n,n'} + \vartheta_t (\rho_m \rho_{m'})^d \log \varrho, \quad t \in [0, d). \qquad (70)$$
Furthermore, if $\Gamma$ is homogeneous then $|\mathbf{m}| - |\mathbf{n}| = |\mathbf{m}'| - |\mathbf{n}'|$, $\rho = \rho^{|\mathbf{m}|-|\mathbf{n}|}$, and, noting that $\rho^d = 1/M$, (70) can be written as

$$I_{m,m'} = \frac{\rho^{-t(|\mathbf{m}|-|\mathbf{n}|)}}{M^2(|\mathbf{m}|-|\mathbf{n}|)} I_{n,n'} + \vartheta_t \frac{(|\mathbf{m}|-|\mathbf{n}|)}{|\mathbf{m}+|\mathbf{m}'|} \log \rho, \quad t \in [0, d). \quad (71)$$

The methodology of [18, 17] involves attempting to

(i) identify a finite number $n_s$ of “fundamental” singular integrals

$$I_{m_{s,1},m'_{s,1},\dots,I_{m_{s,n_s},m'_{s,n_s}},}$$

with $I_{m_{s,1},m'_{s,1}} = I_{\Gamma,\Gamma}$, such that any other singular integral $I_{m,m'}$, for $\mathbf{m}, \mathbf{m}' \in L_h$, is similar to one of them in the sense of Proposition A.1 (via suitable $T, T'$);

(ii) derive an $n_s$-by-$n_s$ linear system of equations satisfied by these fundamental singular integrals, that can be solved to express them (and hence any similar singular integral) in terms of regular integrals that can be computed numerically (e.g. using the composite barycentre rule).

An algorithm for attempting to achieve (i) and (ii) in the general case is presented in [17, Algorithm 1], and the earlier results in [18, §4.3] can be viewed as a specialisation of this algorithm to the case of a disjoint attractor, for which $n_s = 1$. The basic idea is to start from $I_{\Gamma,\Gamma}$ and combine repeated subdivision of the integration domain with repeated applications of Proposition A.1 to determine when integrals are similar. Whether the algorithm succeeds in achieving (i) and (ii) depends on the example being studied, but in [17, §5] it was shown to succeed for a number of well-known examples of non-disjoint attractors including the Sierpinski triangle, Vicsek fractal, Sierpinski carpet and the Koch snowflake. Rather than repeating the full details of the algorithm here, we instead exemplify the procedure in the case of the Koch curve (see Fig. 1(f)), which was not studied in [17].

## A.2 Example - Koch curve

The Koch curve (see Ex. 2.3, Fig. 1(f) and Fig. 3) $\Gamma \subset \mathbb{R}^2$ is the attractor of a homogeneous IFS with $M = 4$ and $\rho = 1/3$, so that $d = \log 4 / \log 3$ and $\rho^d = 1/4$. The only isometries of $\mathbb{R}^2$ under which $\Gamma$ is invariant are the identity, and reflection in the line $x = 1/2$.

To make the notation more compact, given $\mathbf{m} = (m_1, m_2, \dots, m_\ell)$ and $\mathbf{m}' = (m'_1, m'_2, \dots, m'_\ell)$ we write $\Gamma_{\mathbf{m}}, \Gamma_{\mathbf{m}'}$ and $I_{\mathbf{m},\mathbf{m}'}$ as $\Gamma_{m_1 m_2 \dots m_\ell}, \Gamma_{m'_1 m'_2 \dots m'_\ell}$ and $I_{m_1 m_2 \dots m_\ell, m'_1 m'_2 \dots m'_\ell}$, respectively. For example, we write $I_{(1,3),(2,4)}$ as $I_{13,24}$. This compact notation is unambiguous because $M < 10$, which ensures that each entry in $\mathbf{m}$ and $\mathbf{m}'$ is a single-digit integer.
Using Proposition A.1 (with $\varrho = 1$ and appropriate choices of $T$ and $T'$) one can check that $I_{1,1} = I_{2,2} = I_{3,3} = I_{4,4}$, $I_{1,2} = I_{3,4}$, $I_{1,3} = I_{2,4}$, and $I_{i,j} = I_{j,i}$ for all $i,j \in \{1, \dots, 4\}$. Hence a level 1 decomposition of $I_{\Gamma,\Gamma}$ (illustrated in Fig. 3(a)) gives

$$I_{\Gamma,\Gamma} = \sum_{i=1}^{4} \sum_{j=1}^{4} I_{i,j} = 4I_{1,1} + 4I_{1,2} + 2I_{2,3} + R_{\Gamma,\Gamma}, \quad (72)$$

where $R_{\Gamma,\Gamma}$ is a sum of regular integrals, given by

$$R_{\Gamma,\Gamma} := 4I_{1,3} + 2I_{1,4}.$$

The integral $I_{1,1}$ is singular, but by Proposition A.1 is similar to $I_{\Gamma,\Gamma} = I_{0,0}$, with

$$I_{1,1} = \frac{3^t}{16} I_{\Gamma,\Gamma} - \vartheta_t \frac{\log 3}{16}. \quad (73)$$

The integrals $I_{1,2}$ and $I_{2,3}$ are also singular, since $\Gamma_1$ and $\Gamma_2$ intersect at a point, as do $\Gamma_2$ and $\Gamma_3$. But they are not similar to $I_{\Gamma,\Gamma}$, so we apply a level 2 decomposition (see Fig. 3(b)), writing

$$I_{1,2} = I_{14,21} + R_{1,2} \quad \text{and} \quad I_{2,3} = I_{24,31} + R_{2,3}, \quad (74)$$

where $R_{1,2}$ and $R_{2,3}$ are both sums of regular integrals, given by

$$R_{1,2} = 2I_{11,21} + 2I_{11,22} + 2I_{11,23} + I_{11,24} + 2I_{12,21} + 2I_{12,22} + I_{12,23} + 2I_{13,21} + I_{13,22},$$

$$R_{2,3} = 4I_{21,31} + 2I_{21,32} + 2I_{21,33} + 2I_{21,34} + 2I_{22,31} + 2I_{23,31} + I_{23,32}.$$

The integrals $I_{14,21}$ and $I_{24,31}$ are both singular, but are similar to $I_{1,2}$ and $I_{2,3}$, with

$$I_{14,21} = \frac{3^t}{16} I_{1,2} - \vartheta_t \frac{\log 3}{256}, \quad \text{and} \quad I_{24,31} = \frac{3^t}{16} I_{2,3} - \vartheta_t \frac{\log 3}{256}. \quad (75)$$

Combining (72)-(75), we find that the vector of fundamental singular integrals $(I_{\Gamma,\Gamma}, I_{1,2}, I_{2,3})^T$ satisfies the linear system

$$\left( \begin{array}{ccc} \sigma_1 & -4 & -2 \\ 0 & \sigma_2 & 0 \\ 0 & 0 & \sigma_2 \end{array} \right) \left( \begin{array}{c} I_{\Gamma,\Gamma} \\ I_{1,2} \\ I_{2,3} \end{array} \right) = \left( \begin{array}{c} R_{\Gamma,\Gamma} - \vartheta_t \frac{\log 3}{4} \\ R_{1,2} - \vartheta_t \frac{\log 3}{256} \\ R_{2,3} - \vartheta_t \frac{\log 3}{256} \end{array} \right), \quad t \in [0, d), \quad (76)$$

where

$$\sigma_1 = 1 - \frac{3^t}{4}, \quad \sigma_2 = 1 - \frac{3^t}{16}.$$
Solving the system gives

$$I_{2,3} = \frac{1}{\sigma_2} \left( R_{2,3} - \vartheta_t \frac{\log 3}{256} \right), \quad I_{1,2} = \frac{1}{\sigma_2} \left( R_{1,2} - \vartheta_t \frac{\log 3}{256} \right), \quad (77)$$

and

$$I_{\Gamma,\Gamma} = \frac{1}{\sigma_1} \left( R_{\Gamma,\Gamma} + \frac{2}{\sigma_2} \left( 2R_{1,2} + R_{2,3} \right) - \vartheta_t \left( 32 + \frac{3}{\sigma_2} \right) \frac{\log 3}{128} \right). \quad (78)$$

Having derived the representation formulas (77)-(78), one can obtain numerical approximations of $I_{\Gamma,\Gamma}$, $I_{1,2}$ and $I_{2,3}$ by combining these with numerical evaluations of $R_{\Gamma,\Gamma}$, $R_{1,2}$ and $R_{2,3}$, e.g. using the composite barycentre rule with some maximum mesh width $\tilde{h} > 0$.

To apply these results in the assembly of the Galerkin matrix considered in §4(b), we note that since we are using the approximation space $V_N = Y_h$ on the mesh $L_h$ of $\Gamma$, any singular instance of (4.17) will be similar to one of $I_{\Gamma,\Gamma}$, $I_{1,2}$ and $I_{2,3}$, say $I_{i,j}$. To evaluate (4.17) we apply Proposition A.1 to obtain a formula for (4.17) in terms of $I_{i,j}$, which can be evaluated using the value for $I_{i,j}$ already computed (as discussed above). For instance, referring back to Fig. 3(b) with $m = (1, 1, 3)$ and $m' = (1, 1, 4)$ the integral (4.17) is similar to $I_{1,2}$, with

$$I_{113,114} = \frac{3^{2t}}{256} I_{1,2} - \vartheta_t \frac{\log 3}{2048}.$$

To ensure parity between the quadrature accuracies for such singular instances of (4.17), and the regular instances of (4.17), which were computed using the composite barycentre rule with maximum mesh width $h_Q$, we take $\tilde{h} = h_Q/h$.

## Acknowledgements

SC-W was supported by EPSRC grant EP/V007866/1, DH and AG by EPSRC grants EP/S01375X/1 and EP/V053868/1, AM by the PRIN project “NA-FROM-PDEs” and by PNRR-M4C2-I1.4-NC-HPC-Spoke6, funded by the European Union - Next Generation EU, and AC by CIDMA (Center for Research and Development in Mathematics and Applications) and FCT (Foundation for Science and Technology) within project UIDB/04106/2020 (doi.org/10.54499/UIDB/04106/2020). AG, SC-W, DH and AM thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme “Mathematical Theory and Applications of Multiple Wave Scattering”, supported by EPSRC grant EP/R014604/1. AG acknowledges use of the UCL Myriad High Performance Computing Facility (Myriad@UCL) and associated support services. We thank the reviewers for their many helpful comments.
References

[1] NIST Digital Library of Mathematical Functions. dlmf.nist.gov, r1.1.12 of 2023-12-15.

[2] J. BANNISTER, A. GIBBS, AND D. P. HEWETT, *Acoustic scattering by impedance screens/cracks with fractal boundary: well-posedness analysis and boundary element approximation*, Math. Mod. Meth. Appl. Sci. (M3AS), 32 (2022), pp. 291–319.

[3] A. CAETANO, D. P. HEWETT, AND A. MOIOLA, *Density results for Sobolev, Besov and Triebel-Lizorkin spaces on rough sets*, J. Funct. Anal., 281 (2021), p. 109019.

[4] A. M. CAETANO, S. N. CHANDLER-WILDE, A. GIBBS, AND D. P. HEWETT, *Properties of IFS attractors with non-empty interiors and associated function spaces and scattering problems*, In preparation.

[5] A. M. CAETANO, S. N. CHANDLER-WILDE, A. GIBBS, D. P. HEWETT, AND A. MOIOLA, *A Hausdorff-measure boundary element method for acoustic scattering by fractal screens*, Numer. Math., doi.org/10.1007/s00211-024-01399-7, (2024).

[6] M. CEFALO AND M. R. LANCIA, *An optimal mesh generation algorithm for domains with koch type boundaries*, Math. Comput. Simulat., 106 (2014), pp. 133–162.

[7] S. N. CHANDLER-WILDE, I. G. GRAHAM, S. LANGDON, AND E. A. SPENCE, *Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering*, Acta Numer., 21 (2012), pp. 89–305.

[8] S. N. CHANDLER-WILDE AND D. P. HEWETT, *Wavenumber-explicit continuity and coercivity estimates in acoustic scattering by planar screens*, Integr. Equat. Oper. Th., 82 (2015), pp. 423–449.

[9] ——, *Well-posed PDE and integral equation formulations for scattering by fractal screens*, SIAM J. Math. Anal., 50 (2018), pp. 677–717.

[10] S. N. CHANDLER-WILDE, D. P. HEWETT, AND A. MOIOLA, *Sobolev spaces on non-Lipschitz subsets of $\mathbb{R}^n$ with application to boundary integral equations on fractal screens*, Integr. Equat. Operat. Th., 87 (2017), pp. 179–224.

[11] S. N. CHANDLER-WILDE, D. P. HEWETT, A. MOIOLA, AND J. BESSON, *Boundary element methods for acoustic scattering by fractal screens*, Numer. Math., 147 (2021), pp. 785–837.

[12] S. N. CHANDLER-WILDE AND P. MONK, *Wave-number-explicit bounds in time-harmonic scattering*, SIAM J. Math. Anal., 39 (2008), pp. 1428–1455.

[13] X. CLAEYS AND R. HIPTMAIR, *Integral equations on multi-screens*, Integr. Equat. Oper. Th., 77 (2013), pp. 167–197.

[14] L. C. EVANS, *Partial Differential Equations*, AMS, 2010.

[15] L. C. EVANS AND R. E. GARIEPY, *Measure Theory and Fine Properties of Functions*, CRC Press, 2015.
[16] K. FALCONER, *Fractal Geometry: Mathematical Foundations and Applications*, Wiley, 3rd ed., 2014.

[17] A. GIBBS, D. P. HEWETT, AND B. MAJOR, *Numerical evaluation of singular integrals on non-disjoint self-similar fractal sets*, Numer. Alg., doi.org/10.1007/s11075-023-01705-8, (2024).

[18] A. GIBBS, D. P. HEWETT, AND A. MOIOLA, *Numerical evaluation of singular integrals on fractal sets*, Numer. Alg., 92 (2023), pp. 2071–2124.

[19] T. HA-DUONG, *On the transient acoustic scattering by a flat object*, Japan J. Indust. Appl. Math., 7 (1990), pp. 489–513.

[20] ——, *On the boundary integral equations for the crack opening displacement of flat cracks*, Integr. Equat. Oper. Th., 15 (1992), pp. 427–453.

[21] D. P. HEWETT AND A. MOIOLA, *On the maximal Sobolev regularity of distributions supported by subsets of Euclidean space*, Anal. Appl., 15 (2017), pp. 731–770.

[22] P. JONES, J. MA, AND V. ROKHLIN, *A fast direct algorithm for the solution of the Laplace equation on regions with fractal boundaries*, J. Comput. Phys., 113 (1994), pp. 35–51.

[23] A. JONSSON, *Wavelets on fractals and Besov spaces*, J. Fourier Anal. Appl., 4 (1998), pp. 329–340.

[24] A. JONSSON AND H. WALLIN, *Function Spaces on Subsets of $\mathbb{R}^n$, Harwood Academic Publishers, 1984.*

[25] B. N. KHOROMSKIJ AND J. M. MELENK, *Boundary concentrated finite element methods*, SIAM J. Numer. Anal., 41 (2003), pp. 1–36.

[26] J. KRÁL, *Integral operators in potential theory*, Springer, 1980.

[27] P. MATTILA, *Geometry of Sets and Measures in Euclidean Spaces: Fractals and Rectifiability*, CUP, 1995.

[28] V. G. MAZ'YA, *Sobolev Spaces with Applications to Elliptic Partial Differential Equations*, Springer, 2nd ed., 2011.

[29] W. MCLEAN, *Strongly Elliptic Systems and Boundary Integral Equations*, CUP, 2000.

[30] M. MITREA AND M. TAYLOR, *Boundary layer methods for Lipschitz domains in Riemannian manifolds*, J. Funct. Anal., 163 (1999), pp. 181–251.

[31] P. PANAGIOTOPOULOS AND O. PANAGOULI, *The FEM and BEM for fractal boundaries and interfaces. Applications to unilateral problems*, Comput. Struct., 64 (1997), pp. 329–339.

[32] L. RIDDLE, *classic iterated function systems*. larryrid-dle.agnesscott.org/ifs/ifs.htm, downloaded 18 July 2023.

[33] S. A. SAUTER AND C. SCHWAB, *Boundary Element Methods*, Springer, 2011.

[34] H. A. SCHENCK, *Improved integral formulation for acoustic radiation problems*, J. Acoust. Soc. Am., 44 (1968), pp. 41–58.
[35] I. J. ŠNEIBERG, *Spectral properties of linear operators in interpolation families of Banach spaces*, Mat. Issled, 9 (1974), pp. 214–229.

[36] H. TRIEBEL, *Fractals and Spectra*, Birkhäuser, 1997.

[37] T. W. WU AND A. F. SEYBERT, *A weighted residual formulation for the CHIEF method in acoustics*, J. Acoust. Soc. Am., 90 (1991), pp. 1608–1614.