A Hausdorff-measure boundary element method for acoustic scattering by fractal screens

Sound-soft fractal screens can scatter acoustic waves even when they have zero surface measure. To solve such scattering problems we make what appears to be the first application of the boundary element method (BEM) where each BEM basis function is supported in a fractal set, and the integration involved in the formation of the BEM matrix is with respect to a non-integer order Hausdorff measure rather than the usual (Lebesgue) surface measure. Using recent results on function spaces on fractals, we prove convergence of the Galerkin formulation of this ``Hausdorff BEM'' for acoustic scattering in $\mathbb{R}^{n+1}$ ($n=1,2$) when the scatterer, assumed to be a compact subset of $\mathbb{R}^n\times\{0\}$, is a $d$-set for some $d\in (n-1,n]$, so that, in particular, the scatterer has Hausdorff dimension $d$. For a class of fractals that are attractors of iterated function systems, we prove convergence rates for the Hausdorff BEM and superconvergence for smooth antilinear functionals, under certain natural regularity assumptions on the solution of the underlying boundary integral equation. We also propose numerical quadrature routines for the implementation of our Hausdorff BEM, along with a fully discrete convergence analysis, via numerical (Hausdorff measure) integration estimates and inverse estimates on fractals, estimating the discrete condition numbers. Finally, we show numerical experiments that support the sharpness of our theoretical results, and our solution regularity assumptions, including results for scattering in $\mathbb{R}^2$ by Cantor sets, and in $\mathbb{R}^3$ by Cantor dusts.


Introduction
A classical problem in the study of acoustic, electromagnetic and elastic wave propagation is the scattering of a time-harmonic incident wave by an infinitesimally thin screen (or "crack"). In the simplest configuration the incident wave propagates in R n+1 (typically n = 1, 2) and the screen Γ is assumed to be a bounded subset of the hyperplane Γ ∞ = R n × {0}. In standard analyses the set Γ is assumed (either explicitly or implicitly) to be a relatively open subset of Γ ∞ with smooth Figure 1: The first five standard prefractal approximations, Γ 0 , . . . , Γ 4 , to the middle-third Cantor dust Γ, defined by Γ 0 := [0, 1] 2 , Γ n := s(Γ n−1 ), n ∈ N, where s is defined by (4) and (127) with M = 4 and ρ = 1/3. relative boundary ∂Γ. But in a recent series of papers [2,12,14,16] it has been shown how well-posed boundary value problems (BVPs) and associated boundary integral equations (BIEs) for the acoustic version of this screen problem (with either Dirichlet, Neumann or impedance boundary conditions) can be formulated, analysed and discretized for arbitrary screens with no regularity assumption on Γ. In particular, this encompasses situations where either ∂Γ or Γ itself has a fractal nature. The study of wave scattering by such fractal structures is not only interesting from a mathematical point of view, but is also relevant for numerous applications including the scattering of electromagnetic waves by complex ice crystal aggregates in weather and climate science [39] and the modelling of fractal antennas in electrical engineering [45].
Our focus in this paper is on the Dirichlet (sound soft) acoustic scattering problem in the case where Γ itself is fractal. 1 We shall assume throughout that Γ is a compact d-set (see §2.1 for the definition) for some n − 1 < d ≤ n, which in particular implies that Γ has Hausdorff dimension equal to d. More specifically, our attention will be on the special case where Γ is the self-similar attractor of an iterated function system of contracting similarities, in particular on the case where Γ satisfies a certain disjointness condition (described in §2. 3), in which case Γ has (as a subset of R n ) empty interior and zero Lebesgue measure. An example in the case n = 1 is the middle-third Cantor set, which is a d-set for d = log 2/ log 3; an example in the case n = 2 is the middle-third Cantor dust shown in Figure 1, which is a d-set for d = log 4/ log 3. For such Γ, well-posed BVP and BIE formulations for the Dirichlet scattering problem were analysed in [12], where it was shown that the exact solution of the BIE lies in the function space H −1/2 Γ = {u ∈ H −1/2 (Γ ∞ ) : supp u ⊂ Γ} [12, §3.3]. The assumption that d > n −1 implies that this space is non-trivial, and that for non-zero incident data the BIE solution is non-zero, so the screen produces a non-zero scattered field. Our aim in this paper is to develop and analyse a boundary element method (BEM) that can efficiently compute this BIE solution.
One obvious approach, adopted in [16] (and see also [2,29,34]), is to apply a conventional BEM on a sequence of smoother (e.g. Lipschitz) "prefractal" approximations to the underlying fractal screen, such as those illustrated in Figure 1 for the middle-third Cantor dust. (In this example each prefractal is a union of squares.) When Γ has empty interior (as in the current paper) this is necessarily a "non-conforming" approach, in the sense that the resulting discrete approximations do not lie in H −1/2 Γ , the space in which the continuous variational problem is posed. This is because conventional BEM basis functions are elements of L 2 (Γ ∞ ), the intersection of which with H −1/2 Γ is trivial. This complicates the analysis of Galerkin implementations, since Céa's lemma (e.g., [40,Theorem 8.1]), or its standard modifications, cannot be invoked. In [16] we showed how this can be overcome using the framework of Mosco convergence, proving that, in the case of piecewise-constant basis functions, the BEM approximations on the prefractals converge to the exact BIE solution on Γ as the prefractal level tends to infinity, provided that the prefractals satisfy a certain geometric constraint and the corresponding mesh widths tend to zero at an appropriate rate [16,Thm. 5.3]. However, while [16] provides, to the best of our knowledge, the first proof of convergence for a numerical method for scattering by fractals, we were unable in [16] to prove any rates of convergence.
In the current paper we present an alternative approach, in which the fractal nature of the scatterer is explicitly built into the numerical discretization. Specifically, we propose and analyse a "Hausdorff BEM", which is a Galerkin implementation of an H −1/2 Γ -conforming discretization in which the basis functions are the product of piecewise-constant functions and H d | Γ , the Hausdorff d-measure restricted to Γ. A key advantage of the conforming nature of our approximations is that convergence of our Hausdorff BEM can be proved using Céa's lemma. Furthermore, extensions that we make in §3 of the wavelet decompositions from [30] to negative exponent spaces allow us to obtain error bounds quantifying the convergence rate of our approximations, under appropriate and natural smoothness assumptions on the exact BIE solution. While these smoothness assumptions have not been proved for the full range that we envisage (see Proposition 4.9), the convergence rates observed in our numerical results in §6 support a conjecture (Conjecture 4.8) that they hold. Implementation of our Hausdorff BEM requires the calculation of the entries of the Galerkin linear system, which involve both single and double integrals with respect to the Hausdorff measure H d . To evaluate such integrals we apply the quadrature rules proposed and analysed in [25], in which the selfsimilarity of Γ is exploited to reduce the requisite singular integrals to regular integrals, which can be treated using a simple midpoint-type rule. By combining the quadrature error analysis provided in [25] with novel inverse inequalities on fractal sets (proved in §5. 3) we are able to present a fully discrete analysis of our Hausdorff BEM, subject to the aforementioned smoothness assumptions.
An outline of the paper is as follows. In §2 we collect some basic results that will be used throughout the paper on Hausdorff measure and dimension, singular integrals on d-sets, iterated function systems, and function spaces; in particular, in §2. 4 we introduce the function spaces H t (Γ) that are trace spaces on d-sets that will play a major role in our analysis, and recall connections to the classical Sobolev spaces H s Γ established recently in [8]. In §3 we recall from [30] the construction, for n − 1 < d ≤ n, of wavelets on d-sets that are the attractors of iterated function systems satisfying the standard open set condition, and, for n − 1 < d < n, the characterisations of Besov spaces on these d-sets (which we show in Appendix A coincide with our trace spaces H t (Γ) for positive t) in terms of wavelet expansion coefficients. We also extend, in Corollary 3.3, these characterisations, which are crucial to our later best-approximation error estimates, to H t (Γ) for a range of negative t via duality arguments.
In §4 we state the BVP and BIE for the Dirichlet screen scattering problem, showing, in the case when Γ is a d-set, that the BIE can be formulated in terms of a version S of the single-layer potential operator which we show, in Propositions 4.7 and 4.9, maps H t−t d (Γ) to H t+t d (Γ), for |t| < t d and a particular d-dependent t d ∈ (0, 1/2], indeed is invertible between these spaces for |t| < and some 0 < ≤ t d . (The spaces H t d (Γ) ⊂ L 2 (Γ) ⊂ H −t d (Γ) form a Gelfand triple, with L 2 (Γ) the space of square-integrable functions on Γ with respect to d-dimensional Hausdorff measure as the pivot space, analogous to the usual Gelfand triple H 1/2 (Γ) ⊂ L 2 (Γ) ⊂ H −1/2 (Γ) in scattering by a classical screen Γ that is a bounded relatively open subset of Γ ∞ .) Moreover, as Theorem 4.6, we show the key result that, when acting on L ∞ (Γ) (which contains our BEM approximation spaces), S has the usual representation as an integral operator with the Helmholtz fundamental solution as kernel, but now integrating with respect to d-dimensional Hausdorff measure.
In §5 we describe the design and implementation of our Hausdorff BEM, and state and prove our convergence results, showing that, at least in the case that Γ is the disjoint attractor of an iterated function system with n − 1 < dim H (Γ) < n, all the results that are achievable for classical Galerkin BEM (convergence and superconvergence results in scales of Sobolev spaces, inverse and condition number estimates, fully discrete error estimates 2 ) can be carried over to this Hausdorff measure setting (we defer to Appendix B the details of our strongest inverse estimates, derived via a novel extension of bubble-function type arguments to cases where the elements have no interior).
Finally, in §6, we present numerical results, for cases where Γ is a Cantor set or Cantor dust, illustrating the sharpness of our theoretical predictions. We show that our error estimates appear to apply also in cases, such as the Sierpinski triangle, where Γ is not disjoint so that the conditions of our theory are not fully satisfied. We also make comparisons, in terms of accuracy as a function of numbers of degrees of freedom, with numerical results obtained by applying conventional BEM on a sequence of prefractal approximations to Γ, for which we have, as discussed above, only a much more limited theory [16].

Preliminaries
In this section we collect a number of preliminary results that will underpin our analysis.

Hausdorff measure and dimension
For E ⊂ R n and α ≥ 0 we recall (e.g., from [24]) the definition of the Hausdorff α-measure of E, where, for a given δ > 0, the infimum is over all countable covers of E by a collection {U i } i∈N of subsets of R n with diam(U i ) ≤ δ for each i. Where R + := [0, ∞), the Hausdorff dimension of E is then defined to be dim H (E) := sup{α ∈ R + : H α (E) = ∞} = inf{α ∈ R + : H α (E) = 0} ∈ [0, n].
In particular, if E ⊂ R n is Lebesgue measurable then H n (E) = c n |E|, for some constant c n > 0 dependent only on n, where |E| denotes the (n-dimensional) Lebesgue measure of E. Thus dim H (E) = n if E ⊂ R n has positive Lebesgue measure.
As in [31, §1.1] and [43, §3], given 0 < d ≤ n, a closed set Γ ⊂ R n is said to be a d-set if there exist c 2 > c 1 > 0 such that where B r (x) ⊂ R n denotes the closed ball of radius r centred on x. Condition (1) implies that Γ is uniformly locally d-dimensional in the sense that dim H (Γ ∩ B r (x)) = d for every x ∈ Γ and r > 0. In particular (see the discussion in [24, §2.4]) (1) implies that 0 < H d (Γ ∩ B R (0)) < ∞ for all sufficiently large R > 0, so that dim H (Γ) = d.
In the case d = n an open set can satisfy (1). We call Ω ⊂ R n an open n-set if Ω is open and there exists c 1 > 0 such that If Ω is an open n-set then (see [31, Prop. 1 on p. 205]) Γ := Ω is an n-set and |∂Ω| = 0, so that also |∂Γ| = 0. Conversely, if Γ is an n-set and |∂Γ| = 0, it is easy to see that Γ = Ω, where Ω := Γ • , the interior of Γ, and that Ω is an open n-set. We note that open n-sets are also known as interior regular domains, see [8,Def. 5.7].

Singular integrals on compact d-sets
Our Hausdorff BEM involves the discretization of a weakly singular integral equation in which integration is carried out with respect to Hausdorff measure. In order to derive the basic integrability results we require, we appeal to the following lemma, which is [9, Lemma 2.13] with the dependence of the equivalence constants made explicit.

Iterated function systems
The particular example of a d-set we focus on in this paper is the attractor of an iterated function system (IFS) of contracting similarities, by which we mean a collection {s 1 , s 2 , . . . , s M }, for some M ≥ 2, where, for each m = 1, . . . , M , s m : R n → R n , with |s m (x) − s m (y)| = ρ m |x − y|, x, y ∈ R n , for some ρ m ∈ (0, 1). The attractor of the IFS is the unique non-empty compact set Γ satisfying We shall assume throughout that Γ satisfies the open set condition (OSC) [24, (9.11)], meaning that there exists a non-empty bounded open set O ⊂ R n such that Then [43,Thm. 4 For a homogeneous IFS, where ρ m = ρ ∈ (0, 1) for m = 1, . . . , M , the solution of (6) is Returning to the general, not necessarily homogeneous, case, the OSC (5) also implies (again, see [43,Thm. 4.7]) that Γ is self-similar in the sense that the sets which are congruent copies of Γ, satisfy That is, Γ can be decomposed into M congruent copies of itself, whose pairwise intersections have Hausdorff measure zero. For many of our results we make the additional assumption that the sets Γ 1 , . . . , Γ M are disjoint, in which case we say that the IFS attractor Γ is disjoint. We recall that if Γ is disjoint then it is totally disconnected [24,Thm. 9.7]. Examples of IFS attractors that are disjoint are the Cantor set and Cantor dust (Figure 1), while an example that satisfies the OSC but is not disjoint is the Sierpinski triangle. The next result relates disjointness to the OSC.
The following lemma, which shows that Γ is disjoint only if d < n, motivates the restriction of our results in large parts of §3 to the case d < n.

Function spaces on subsets of R n
Here we collect some results on function spaces from [8,14,30]. Our function spaces will be complexvalued, and we shall repeatedly use the following terminology relating to dual spaces 3 . If X and Y are Hilbert spaces, X * is the dual space of X, and I : Y → X * is a unitary isomorphism, we say that (Y, I) is a unitary realisation of X * and define the duality pairing y, x Y ×X := Iy(x), for y ∈ Y, x ∈ X. Having selected a unitary realisation (Y, I) of X * , we adopt (X, I * ) as our unitary realisation of Y * , where I * x(y) := y, x Y ×X , so that x, y X×Y = y, x Y ×X , for y ∈ Y, x ∈ X.
For s ∈ R let H s (R n ) denote the Sobolev space of tempered distributions ϕ for which the norm 4 . We recall that (H s (R n )) * can be unitarily re- and ψ ∈ H s (R n ), so that the resulting duality pairing ·, · H −s (R n )×H s (R n ) extends both the L 2 (R n ) inner product and the action of tempered distributions on Schwartz functions (see, e.g., [14, §3.1.3]). For an open set Ω ⊂ R n we denote by H s (Ω) the closure of C ∞ 0 (Ω) in H s (R n ), and for a closed set E ⊂ R n we denote by H s E the set of all elements of H s (R n ) whose distributional support is contained in E. These two types of spaces are related by duality, with (H −s E , I) being a unitary realisation of ( . Note also that, if E ⊂ Ω ⊂ R n and E is compact, Ω is open, then, as a consequence of [32,Lemma 3.24], H s E is a closed subspace of H s (Ω). In addition to the spaces just introduced we use, at some points, the standard Sobolev space H s (Ω), for Ω ⊂ R n open and s ∈ R, defined as the space of restrictions to Ω of the distributions ϕ ∈ H s (R n ), equipped with the quotient norm u H s (Ω) := inf ϕ∈H s (R n ) ϕ| Ω =u ϕ H s (R n ) (see [32], [14, §3.1.4]).
Fix 0 < d ≤ n and let Γ ⊂ R n be H d -measurable. We denote by L 2 (Γ) the space of (equivalence classes of) complex-valued functions on R n that are measurable and square integrable with respect to H d | Γ , normed by In practice we shall view L 2 (Γ) as a space of functions on Γ, by identifying elements of L 2 (Γ) with their restrictions to Γ. The dual space (L 2 (Γ)) * can be realised in the standard way as (L 2 (Γ), I), where I : . Now assume that Γ is a d-set and that 0 < d ≤ n. Then function spaces on R n and Γ are related via the trace operator tr Γ of [43, §18.5]. Defining tr Γ (ϕ) = ϕ| Γ ∈ L 2 (Γ) for ϕ ∈ C ∞ 0 (R n ), one can show [43,Thm 18.6 which we assume through the rest of this section, and 0 < d < n, then tr Γ extends to a continuous linear operator tr Γ : H s (R n ) → L 2 (Γ) with dense range. 5 This trivially holds also for d = n, since the embedding of H s (R n ) into L 2 (R n ), for s > 0, and the trace tr Γ : L 2 (R n ) → L 2 (Γ) are both continuous with dense range. Setting we define the trace space H t (Γ) := tr Γ (H s (R n )) ⊂ L 2 (Γ), which we equip with the quotient norm This makes H t (Γ) a Hilbert space unitarily isomorphic to the quotient space H s (R n )/ ker(tr Γ ). Clearly, for t > t > 0, and the embeddings are continuous with dense range. As explained in [8,Rem. 6.4] [8, §6]), for the case 0 < d < n, under the further assumption that t < 1, H t (Γ) coincides (with equivalent norms) with the Besov space B t 2,2 (Γ) of [31]. Arguing in the same way, using [31, Theorem VI.1] and that H s (R n ) coincides with the Besov space B s 2,2 (R n ) (e.g., [31, p. 8]) 6 , this holds also for d = n.
In the particular case that Γ = Ω and Ω is an open n-set in R n , it holds that 7 in the sense that the operation of restriction to Ω establishes a one-to-one correspondence from H s (Γ) onto H s (Ω), the norms of corresponding elements being the same.
For t > 0 we denote by H −t (Γ) the dual space (H t (Γ)) * . Since (11) holds, and the embeddings are continuous with dense range, also H −t (Γ) ⊂ H −t (Γ), for t > t > 0, and this embedding is continuous with dense range. Further, via the Riesz map I : L 2 (Γ) → L 2 (Γ) * introduced above, L 2 (Γ) is continuously and densely embedded in H −t (Γ), for t > 0. Setting H 0 (Γ) = L 2 (Γ), and combining these embeddings, we then have that H t (Γ) is embedded in H t (Γ) with dense image for any t, t ∈ R with t > t, and that if g ∈ H t (Γ) for some t ≥ 0 and f ∈ L 2 (Γ) then Suppose that (10) holds. By the definition of H t (Γ), tr Γ : H s (R n ) → H t (Γ) is a continuous linear surjection with unit norm. Its Banach space adjoint tr Γ * : is then a continuous linear injection with unit norm, and composing tr Γ * with the unitary isomorphism (I −s ) −1 produces a continuous linear in- For 0 < d < n the operator trΓ extends further to a surjective map trΓ : B is dense), but we do not need that generality here. For details see [8, §6] and [43,Thm 18.6]. 6 Our standard notation for Besov spaces on R n is that of, e.g., [44]. 7 To see this, note that H s (Γ) and H s (Γ) are unitarily isomorphic to the quotient spaces H s (R n )/ ker(trΓ) and H s (R n )/ ker(|Ω), respectively, where trΓ : H s (R n ) → L2(Γ) and |Ω : H s (R n ) → L2(Ω), and that, if |∂Ω| = 0 (which holds if Ω is an open n-set), then these quotient spaces coincide (with equality of norms) since ker(trΓ) = ker(|Ω). This last claim follows from the fact that the operation of restriction to Ω is an isomorphism from L2(Γ) onto L2(Ω) if |∂Ω| = 0 (recall that H n is a constant times the Lebesgue measure). Figure 2: The main function spaces introduced in §2.4 and their relations. Here Γ ⊂ R n is a d-set with 0 < d ≤ n, t = s − n−d 2 ∈ (0, 1), both arrows represent unitary isomorphisms. Where we write A ⊂ B, the function space A is densely and continuously embedded in B. Where ∩ separates A and B vertically, the space A is a closed subspace of B.
with unit norm, which satisfies 8 In particular, when f ∈ L 2 (Γ) we have that Since Lemma 3.2]. Key to our analysis will be the following stronger result. This is proved, for the case 0 < d < n, in [8, Prop. 6.7, Thm 6.13], and the arguments given there (we need only the simplest special case m = 0) extend to the case d = n, with the twist that, to justify the existence of a bounded right inverse E Γ,0 in Step 2 of the proof of [8, Prop. 6.7], we need to use (as above) that H s (R n ) = B s 2,2 (R n ) and [31, Thm. VI.3 on p. 155]. Figure 2 shows the main relations between these function spaces.
is a unitary isomorphism. Accordingly, the range of tr * Γ is equal to H −s Γ , and Furthermore, tr * Γ (L 2 (Γ)) is dense in H −s Γ .

Function spaces on screens
For the screen scattering problem, we define function spaces on the hyperplane Γ ∞ := R n × {0} and subsets of it (for example, on the compact subset Γ ⊂ Γ ∞ that forms the screen) by associating Γ ∞ 8 We omit in our notation for trΓ : H s (R n ) → H t (Γ) any dependence on s. This is justified since, for every s > n−d 2 , trΓϕ = ϕ|Γ for ϕ ∈ C ∞ 0 (R n ), and C ∞ 0 (R n ) is dense in H s (R n ). Likewise, we omit any dependence on t in our notation for tr * Γ : To see that this is justified, in particular that the values of the tr * Γ operators coincide where their domains intersect, denote tr * Γ : H −t (Γ) → H −s (R n ) temporarily by tr * Γ,t to make explicit the domain. Then, for t > t > 0, tr * Γ,t f = tr * Γ,t f , for f ∈ H −t (Γ) ⊂ H −t (Γ), as a consequence of (14).

Wavelet decompositions
The spaces H t (Γ) were defined in §2.4 for Γ ⊂ R n a general d-set with 0 < d ≤ n. In the case where Γ is a disjoint IFS attractor (in the sense of §2.3, in which case, by Lemma 2.6, d < n), it was shown in [30] that, when n − 1 < d = dim H (Γ) < n, the spaces H t (Γ) can be characterised in terms of wavelet decompositions for t > 0. 10 This, more precisely an extension of this characterisation to negative t, will be central to our BEM convergence analysis later. In the current section we recap the notation and main results from [30] that we will need, initially assuming only that n − 1 < d ≤ n and that the OSC is satisfied.
Let Γ be the attractor of an IFS {s 1 , . . . , s M } as in (4), and assume that the OSC (5) holds. Following An example of use of the notation E m is shown in Figure 3. Note that this notation extends that of (8) where the sets Γ 1 , . . . , Γ M were introduced, corresponding to the case E = Γ and = 1 here. 9 Our notation follows [32]. Given an open set Ω ⊂ R n+1 , W 1 (Ω) is the Sobolev space whose norm is defined "intrinsically", via integrals over Ω, while H 1 (Ω) is defined (see §2.4) "extrinsically" as the set of restrictions to Ω of functions in H 1 (R n+1 ). These spaces coincide if Ω is Lipschitz (e.g., [32,Theorem 3.16]), in particular if Ω = R n+1 , but not in general, in particular not, in general, for Ω = D. 10 More precisely, [30] showed that the Besov spaces B p,q α (Γ), as defined in [30, §6] (and see Definition A.1 below), can be characterised in this way for α > 0 and 1 ≤ p, q ≤ ∞, provided Γ preserves Markov's inequality, which holds in particular (see Remark A.6) if d > n − 1. We show in Appendix A (see Corollary A.4 and Remark A.6) that, for t > 0, H t (Γ) = B 2,2 t (Γ) with equivalence of norms if Γ is a d-set with d > n − 1, so that this characterisation carries over to our trace spaces H t (Γ). Let W 0 be the space of constant functions on Γ, a one-dimensional subspace of L 2 (Γ) spanned by More generally, for ∈ N let Since H d (Γ m ∩Γ m ) = 0 for m = m (a consequence of (9) (self-similarity)), W is a M -dimensional subspace of L 2 (Γ) with orthonormal basis Clearly W 0 ⊂ W 1 ⊂ W 2 ⊂ · · · ⊂ L 2 (Γ). But the bases we introduced above are not hierarchical, in the sense that the basis for W +1 does not contain that for W . Following [30] we introduce hierarchical wavelet bases on the W spaces by decomposing where W +1 W denotes the orthogonal complement of W in W +1 . As already noted, W 0 is one-dimensional, with orthonormal basis For ∈ N the space W +1 W is ((M − 1)M )-dimensional, and an orthonormal basis of W +1 W (see Figure 4)   Hence, recalling (19) and setting ψ m 0 = ψ m for m = 1, 2, . . . , M − 1, we obtain the following orthonormal basis of W : These bases are hierarchical, in the sense that the basis for W +1 contains that for W . Furthermore, as noted in [30, p. 334] (and demonstrated in the proof of Theorem 5.1 below), elements of L 2 (Γ) can be approximated arbitrarily well by elements of W as → ∞, which implies that is a complete orthonormal set in L 2 (Γ). Hence every f ∈ L 2 (Γ) has a unique representation with and The next result, which combines Theorems 1 and 2 in [30], provides a characterization of the space H t (Γ) ⊂ L 2 (Γ) (introduced after (10)) for n − 1 < d = dim H (Γ) < n and 0 < t < 1 in terms of the wavelet basis introduced above, under the assumption that Γ is a disjoint IFS attractor. This assumption ensures that the piecewise-constant spaces W are contained in H t (Γ) for all t > 0 (see [30,Thm. 2]). In the statement of Theorem 3.1 the set J ν is defined for ν ∈ Z by and ν 0 is defined to be the unique integer such that 0 ∈ J ν 0 , i.e. such that 2 −ν 0 ≤ diam(Γ) < 2 −ν 0 +1 . Note that 11 with disjoint unions on both sides, so m∈I F (m) whenever the convergence is unconditional, as is the case, for instance, for (21). For convenience we introduce in this theorem a norm · t that is different, but trivially equivalent to that used in [30], which was ). Let Γ be a disjoint IFS attractor with n − 1 < d = dim H (Γ) < n, and let 0 < t < 1.
Theorem 3.1 has the following important corollary, which is obtained by duality. In this corollary and subsequently (see, e.g., [13,Remark 3.8]), given an interval I ⊂ R we will say that a collection of Hilbert spaces {H s : s ∈ I}, indexed by I, is an interpolation scale if, for all s, t ∈ I and 0 < η < 1, with equivalent norms. We will say that {H s : s ∈ I} is an exact interpolation scale if, moreover, the norms of (H s , H t ) η and H θ coincide, for all s, t ∈ I and 0 < η < 1. 11 As an example of these definitions, suppose that Γ is the (disjoint) attractor of the IFS illustrated in Figure 3, in which case diam(Γm) = diam(Em), where Em is as defined in Figure 3, in particular diam(Γ0) = diam(Γ) = 1. Then 12 Here, and subsequently, (Hs, Ht)η denotes the standard complex interpolation space, (Hs, Ht) [η] in the notation of [3]; equivalently, the K-or J-method real interpolation spaces denoted (Hs, Ht)η,2 in [3,13], which are the same Hilbert spaces (Hs, Ht) [η] , with equal norms, if the K-and J-methods are appropriately normalised (see [13,Remark 3.6] and [15]). Corollary 3.3. Let Γ and t satisfy the assumptions of Theorem 3.1.
(iv) The duality pairing ·, · H −t (Γ)×H t (Γ) can be evaluated using the wavelet basis as . With this pairing, H −t (Γ) provides a unitary realisation of (H t (Γ)) * with respect to the norms · t and · −t .
Proof. For t ∈ R we can define the weighted 2 sequence space which is a Hilbert space with the obvious inner product. The dual space of h t can be unitarily realised as h −t , with duality pairing Furthermore, Theorem 3.1 implies that the space H t (Γ) is linearly and topologically isomorphic to h t for 0 < t < 1 (unitarily if we equip H t (Γ) with · t ). Hence by duality H −t (Γ) is also linearly and topologically isomorphic to h −t for the same range of t (unitarily if we equip H −t (Γ) with · −t ). From these observations parts (i)-(iv) of the result follow, noting, in the case of (iv), that f is a continuous antilinear functional on H t (Γ).
For (v) we note that by [ . . , M − 1}}, µ to be the counting measure on X , A to be the map taking f ∈ H τ j (Γ) to the sequence β defined by (22) or (26) (as appropriate), and A basic interpolation result is that if X 0 ⊃ X 1 and Y 0 ⊃ Y 1 are Hilbert spaces, with X 1 and Y 1 continuously embedded in X 0 and Y 0 , respectively, and A : X j → Y j is a linear and topological isomorphism, for j = 0, 1, then, for 0 Remark 3.4. Theorem 3.1 and Corollary 3.3, together with the above interpolation result applied with A taken as the identity operator, imply that, if Γ satisfies the assumptions of Theorem 3.1, then

or indeed any other equivalent norm).
We include the following corollary, although we will not use it subsequently, because it may be of independent interest. Note that the range of s does not extend to s = −(n − d)/2 since, by [ for all s ∈ R and 0 < θ < 1. Proof. Apply the basic interpolation result above with A = tr * Γ and , where s and t are related by (16) and similarly t = s − (n − d)/2; note that A : X j → Y j is a linear and topological isomorphism, for j = 0, 1, by Theorem 2.7. This gives, for 0

BVPs and BIEs
In this section we state the BVP and BIE that we wish to solve. We consider time-harmonic acoustic scattering of an incident wave u i propagating in R n+1 (n = 1, 2) by a planar screen Γ, a subset of the hyperplane Γ ∞ = R n × {0}. We initially consider the case where Γ is assumed simply to be non-empty and compact, for which a well-posed BVP/BIE formulation was presented in [16, §3.2].
We later specialise to the case where Γ is a d-set for some n − 1 < d ≤ n, and then further to the case where Γ is a disjoint IFS attractor.
Our BVP, stated as Problem 4.1 below, is for the scattered field u, which is assumed to satisfy the Helmholtz equation in D := R n+1 \ Γ, for some wavenumber k > 0, and the Sommerfeld radiation condition We assume that the incident wave u i is an element of W 1,loc (R n+1 ) satisfying (27) in some neighbourhood of Γ (and hence C ∞ in that neighbourhood by elliptic regularity, see, e.g., [22,Thm 6 , which motivates the following definition, in which P denotes the orthogonal projection Problem 4.1. Let Γ ⊂ Γ ∞ be non-empty and compact. Given k > 0 and g ∈ H 1/2 (Γ c ) ⊥ , find u ∈ C 2 (D) ∩ W 1,loc (D) satisfying (27) in D, (28), and the boundary condition for some (and hence every) σ ∈ C ∞ 0,Γ . In the case of scattering of an incident wave u i , g is given specifically as The next result reformulates the BVP as a BIE. In this theorem S : H where φ denotes the complex conjugate of φ, 15 is the Hankel function of the first kind of order zero (e.g., [1,14 Here and in what follows, Γ c will denote Γ∞ \ Γ, the complement of Γ in Γ∞ (rather than its complement in R n+1 , which we have denoted by D). 15 If φ ∈ L2(Γ∞) ⊂ H −1/2 (Γ∞),φ is the usual complex conjugate, and this definition of the complex conjugate is extended to H −1/2 (Γ∞) by density. Equation (9.1.3)]), and σ is any element of C ∞ 0,Γ with x ∈ supp σ. The ± in (32) indicates that either trace can be taken, with the same result. In the case when Γ is the closure of a Lipschitz open subset of Γ ∞ and φ ∈ L 2 (Γ) the potential can be expressed as an integral with respect to (Lebesgue) surface measure, namely The operator S : H where σ ∈ C ∞ 0,Γ is arbitrary, which is continuous and coercive (see Lemma 4.3 below).
with g given by (31) in the case of scattering of an incident wave u i .
Define the sesquilinear form a(·, ·) on H Then the BIE (36) can be written equivalently in variational form as: given This equation will be the starting point for our Galerkin discretisation in §5.
The following lemma provides conditions under which a compact screen Γ ⊂ Γ ∞ produces a nonzero scattered field. We note that a sufficient condition for H Thm 2.12], and that when Γ is a d-set this is also a necessary condition [27,Thm 2.17].

The BIE on d-sets in trace spaces
Suppose now that Γ is a compact d-set with n − 1 < d ≤ n. The assumption that d > n − 1 ensures that Γ produces a non-zero scattered field under the conditions of Lemma 4.4. Furthermore, the condition d > n − 1 is equivalent to (10) with s = 1/2, so that the results in §2.4 apply with s = 1/2 and It then follows from (32) and (15) that for Ψ ∈ L 2 (Γ) the potential S has the following integral representation with respect to Hausdorff measure: is well-defined for all x ∈ R n+1 and is continuous, and Proof. It is clear that F (x) is well-defined for x ∈ D since the integrand is then in L ∞ (Γ), as Φ(x, y) is continuous for x = y, and H d (Γ) < ∞ as Γ is a bounded d-set. For some constant C > 0, whereê ∈ R n+1 is any unit vector, and let F ε (x) : is continuous, in the (usual) sense that it is equal almost everywhere with respect to n + 1-dimensional Lebesgue measure to a continuous function (the function F ), by (45).
Noting that tr Γ : then S is continuous and coercive, with the same associated constants as S (the constants in Lemma 4.3). Furthermore, it follows from (37) and (14) that and the variational problem (38) can be equivalently stated as: and the solutions of (49) and (38) are related through φ = tr * Γ Ψ. A schematic showing the relationships between the relevant function spaces and operators is given in Figure 5.
The following integral representation for S will be crucial for our Hausdorff BEM in §5. Figure 5: Schema of relevant function spaces and operators for s = 1/2, t = t d : , and the result follows.
The definition and mapping properties of tr Γ and tr * Γ , noted in §2.4, combined with the representation (50), enable us to extend the domain of S to H −t (Γ), for t d < t < 2t d , or restrict it to H −t (Γ), for 0 < t < t d , as stated in the following key lemma.
and is continuous.
The claimed mapping property of S follows from (50) since , with s and t related by (10), and since the mapping We now make a conjecture concerning the mapping properties of S −1 . To the best of our knowledge there are no results in this direction in the case d < n, but the conjecture can be seen as an extension of known results in the case d = n, since the conjecture is known to be true in the case that Γ = Ω for some bounded Lipschitz domain Ω ⊂ Γ ∞ , in which case d = n (so t d = 1/2) and Ω is an open n-set so that (12) applies. For in this case the single layer BIO S Ω , defined below (39), is invertible as an operator from H t−1/2 (Ω) to H t+1/2 (Ω) for |t| < 1/2 -see [42,Thm 1.8] for the case n = 1 and [36, Theorem 4.1] 16 for the case n = 2 -which implies that S : H t−t d (Γ) → H t+t d (Γ) is invertible, by (12), Theorem 2.7, and the fact that H s (Ω) = H s Ω (see, e.g., [32,Thm 3.29]). A further motivation for our conjecture is that its truth implies convergence rates for our BEM (defined and analysed in §5) that are evidenced by numerical experiments in §6 for cases with n − 1 < d < n.
is invertible, and hence (by Proposition 4.7 and the bounded inverse theorem) a linear and topological isomorphism.
While we are not able to prove Conjecture 4.8 in its full generality, in the case where Γ is a disjoint IFS attractor, in which case, by Lemma 2.6, d < n, we can prove that S is invertible for a range of t, using Corollary 3.3 and results from function space interpolation theory.
Proof. We recall the following result from [33,Prop. 4.7], which quotes [38]. Suppose E j and F j are Banach spaces, E 1 ⊂ E 0 and F 1 ⊂ F 0 with continuous embeddings, and T : E j → F j is a bounded linear operator for j = 0, 1. Let E θ = (E 0 , E 1 ) θ be defined by complex interpolation, and similarly F θ , so T : E θ → F θ and is bounded, for 0 < θ < 1. Assume that T : in a neighbourhood of t = 0 then follows by taking 0 < τ < t d and applying the above result with is invertible and tr Γ g ∈ H t+t d (Γ), then, again using the mapping properties of tr * Γ from Theorem 2.7, the solution φ = tr * Γ Ψ of the BIE (36) satisfies for some constant C > 0 independent of φ and g. In the case of scattering of an incident wave u i , in which g is given by (31), we have that tr Γ g = −tr Γ γ ± ((σu i )| U ± ) ∈ H t+t d (Γ) for all 0 < t < t d , since u i is C ∞ in a neighbourhood of Γ. Hence, in this case, if Γ is a disjoint IFS attractor with n − 1 < d = dim H (Γ) < n, then (53) holds for 0 < t < , where is as in Proposition 4.9, and, if Conjecture 4.8 holds, then (53) holds for 0 < t < t d whenever Γ is a compact d-set with n − 1 < d ≤ n.

The Hausdorff BEM
We now define and analyse our Hausdorff BEM. To begin with, we assume simply that Γ is a compact d-set for some n − 1 < d ≤ n.
Given N ∈ N let {T j } N j=1 be a "mesh" of Γ, by which we mean a collection of H d -measurable subsets of Γ (the "elements") such that H d (T j ) > 0 for each j = 1, . . . , N , H d (T j ∩ T j ) = 0 for j = j , and Define the N -dimensional space of piecewise constants and set Our proposed BEM for solving the BIE (36), written in variational form as (38), uses V N as the approximation space in a Galerkin method. Given g ∈ ( H 1/2 (Γ c )) ⊥ we seek φ N ∈ V N such that (with a defined by (37)) be the corresponding basis for V N . Then, writing φ N = N j=1 c j e j , (56) implies that the coefficient vector c = (c 1 , . . . , c N ) T ∈ C N satisfies the system where, by (48), (13), and (51), the matrix A ∈ C N ×N has (i, j)-entry given by and, by (15), the vector b ∈ C N has ith entry given by For the canonical L 2 (Γ)-orthonormal basis for V N , where f j | T j := (H d (T j )) −1/2 and f j | Γ\T j := 0, j = 1, . . . , N , Once we have computed φ N by solving (57) we will compute approximations to u(x) and u ∞ (x), given by (35)/ (32) and (42), respectively. Each expression takes the form 17 for some ϕ ∈ ( H 1/2 (Γ c )) ⊥ . Explicitly, where σ is any element of C ∞ 0,Γ (with x not in the support of σ in the case u(x)), In each case we approximate J(φ) by J(φ N ) which, recalling (15), is given explicitly by For the canonical L 2 (Γ)-orthonormal basis for V N , The following is a basic convergence result.
Theorem 5.1. Let Γ be a compact d-set for some n − 1 < d ≤ n. Then for each N ∈ N the variational problem (56) has a unique solution Proof. The well-posedness of (56) follows from the Lax-Milgram lemma and the continuity and coercivity of S : H −1/2 Γ → ( H 1/2 (Γ c )) ⊥ . Furthermore, by Céa's lemma (e.g., [40,Theorem 8.1]) we have the following quasi-optimality estimate, where C a , α > 0 are as in (41): 17 Note that if ψ ∈ H Hence to prove convergence of φ N to φ it suffices to show that inf ψ N ∈V N φ − ψ N H −1/2 (Γ∞) → 0 as N → ∞, which, by the definition of V N and the fact that tr * Γ : To show the latter we note that the space C(Γ) of continuous functions on Γ (equipped with the supremum norm) is continuously embedded in L 2 (Γ) with dense image (continuity is obvious and density follows by the density of C ∞ 0 (R n ) in H 1/2 (R n ) and that (see §2.4) tr Γ : H 1/2 (R n ) → L 2 (Γ) is continuous and has dense range). Then, given > 0 and Ψ ∈ L 2 (Γ), there exists Ψ ∈ C(Γ) such that Ψ − Ψ L 2 (Γ) < /2, and by the uniform continuity of Ψ and the fact that is a bounded linear functional.

Best approximation error estimates
We now assume that Γ is the attractor of an IFS of contracting similarities satisfying the OSC, as in §2.3. In this case, one possible choice of BEM approximation space could be V N = W , for some ∈ N 0 (recall the definitions in §3), so that {T j } N j=1 = {Γ m } m∈I and N = M . However, since the spaces W are defined by refinement to a certain prefractal level, if the contraction factors ρ 1 , . . . , ρ M are not all equal the resulting mesh elements may differ significantly in size when is large. This motivates the use of spaces defined by refinement to a certain element size. Recalling the notations of §3 (in particular, J ν is defined in (23) and ν 0 below (23)), for ν ≥ ν 0 let where Note that both of these spanning sets for X ν are orthonormal bases, and that X ν ⊂ X ν , for ν 0 ≤ ν ≤ ν , giving that (but with this spanning set not linearly independent). Note also that, for ν ≥ ν 0 , {Γ m : m ∈ K ν } is a mesh in the sense introduced at the beginning of this section. Figure 6 illustrates the definitions of X ν and K ν , showing the meaning of K 2 for a particular IFS for which diam(Γ) = 1 so that ν 0 = 0; for this same IFS we have that K 0 = K 1 = {1, 2, 3}.
Furthermore, for −1 < t 1 < t 2 < 1 there exists c > 0, independent of ν and f , such that Proof. The boundedness result follows from Theorem 3.1 and Corollary 3.3, noting that For the approximation result, we note, where c > 0 denotes a constant independent of ν and f , not necessarily the same at each occurrence, that Let X ν := tr * Γ (X ν ). (Recall that the choice of s > (n − d)/2 in the definition of tr * Γ makes no difference to the definition of X ν .) Since X ν ⊂ L 2 (Γ) we have X ν ⊂ H s Γ for all s < −(n − d)/2. Then the following best approximation error estimate follows immediately, on application of tr * Γ , from Theorem 2.7 and Proposition 5.2 applied with −1 < t 1 < t 2 < 0. Corollary 5.3. Let Γ be a disjoint IFS attractor with n − 1 < d = dim H (Γ) < n, and suppose that −(n − d)/2 − 1 < s 1 < s 2 < −(n − d)/2. Then there exists c > 0, independent of ν and ψ, such that Remark 5.4. The above corollary holds also for larger values of s 2 , but is then a trivial result since, as noted above Corollary 3.5, H s 2 Γ = {0} for s 2 ≥ −(n − d)/2.
Corollary 5.3 can be rephrased with 2 −ν replaced by a more general element size h, and X ν replaced by where and This framework is a natural one for numerics -we specify h > 0 and consider all those components which have diameter less than or equal to h but whose "parent" (in the IFS structure) has diameter greater than h. The change from < to ≤ and ≥ to > in moving from K ν to L h is intentional, so that the definition of L h allows components that are equal to h in diameter; h is an upper bound for, and may be equal to, the maximum element diameter (this spanning set not linearly independent). We assume in the following corollary and subsequently that h, an upper bound for h N ≤ diam(Γ), satisfies Except where indicated explicitly otherwise, results through the rest of the paper hold for all h in this range. The following corollary, a best approximation result for the spaces Y h , follows from the inclusion X ν ⊂ Y h , for some suitable ν, and Corollary 5.3.

Corollary 5.5. Under the assumptions of Corollary 5.3 we have that
for some constant c > 0 independent of h and ψ.
Proof. It is enough to show (79) for 0 < h ≤ diam(Γ)/2, since the result trivially holds for larger h (just take ψ h = 0 on the left hand side). Thus, given 0 so that, where c > 0 denotes a constant independent of h and ψ, not necessarily the same at each occurrence,

Galerkin error estimates
We now use the best approximation error results proved above to give an error bound for the Galerkin approximation to the BIE (36) when the solution φ is sufficiently smooth. Note that, combining the results of this theorem with inverse estimates that we prove as Theorem 5.10 below, we extend (80) to a bound on φ − φ N H s 1 Γ for a range of s 1 in Corollary 5.12 below.
Proof. We first note that, in the general notation introduced at the start of this section, to get Then to obtain (80) we simply combine the best approximation error bound (79) (for s 1 = −1/2 and s 2 = s) with the quasioptimality estimate (66).
Remark 5.8 (Convergence rates when Γ is a disjoint homogeneous IFS). Suppose that, in addition to the assumptions of Theorem 5.6, Γ is homogeneous, with ρ m = ρ for m = 1, . . . , M , for some 0 < ρ < 1. Then, taking h = ρ diam Γ, i.e. using the approximation space V N = tr * Γ (span ({χ m } m∈I )), Theorem 5.6 implies that (Here and below c > 0 denotes some constant independent of , φ, and ζ, not necessarily the same at each occurrence.) If φ ∈ H s Γ for all s < −(n − d)/2 we can take s = −(n − d)/2 − for arbitrarily small in the first of the above estimates, and the same holds for the second of the above estimates if also ζ satisfies ζ ∈ H s Γ for the same range of s. (If Conjecture 4.8 holds this smoothness of φ and ζ is guaranteed for all scattering problems by Remark 5.7 if ϕ satisfies the conditions in that remark.) Then, recalling that d = log(1/M )/ log ρ, we get s + 1/2 = (d + 1 − n)/2 − = (1/2)(log(1/M )/ log ρ + 1 − n − 2 ), so that, for each > 0, In the case n = 1 (e.g. Γ a Cantor set, see (126), for which M = 2), the fact that > 0 can be taken arbitrarily small means that in numerical experiments we expect to see errors in computing φ N and J(φ N ) that tend to zero roughly like M − /2 and M − respectively, independent of the parameter ρ, if the above bounds are sharp.
In the case n = 2 (e.g. Γ a Cantor dust, see (127), for which M = 4), we expect errors in computing φ N and J(φ N ) that tend to zero roughly like (M ρ) − /2 and (M ρ) − respectively, if the above bounds are sharp. Note that in this case we need ρ > 1/M to ensure d = log(1/M )/ log(ρ) > 1 = n − 1, and that these predicted convergence rates, as a function of , decrease as d approaches 1 with M fixed. For the Cantor dust with ρ = 1/3 (the "middle-third" case) we predict convergence rates of roughly (3/4) /2 and (3/4) , respectively.
Remark 5.9 (Connection to standard BEM convergence results). The results of Theorem 5.6, because they require that d < n, do not apply when Γ is the closure of a Lipschitz domain (so that d = dim H (Γ) = n), for which case standard regularity and convergence results (e.g., [21,41,42] and see the discussion above Conjecture 4.8) predict that φ ∈ H s Γ and that the bounds (80) and (81) hold for all s < 0. Note that these convergence rates for standard BEM are those predicted by taking the formal limit as d → n − in the results of Theorem 5.6. In §6 we will compare results for the cases where Γ is a Cantor set or Cantor dust, the attractor of the IFS (126) or (127), respectively, with results for standard BEM. If we take the parameter ρ = 0.5 in each of (126) and (127) then the attractor is just Γ = [0, 1] n , with n = 1 or 2. As we note in §6.1 and §6.2, the relative errors φ − φ N H −1/2 Γ / φ H −1/2 Γ for our new Hausdorff-measure BEM for the Cantor set and Cantor dust with ρ = 0.49 are almost the same as relative errors for the limiting case ρ = 0.5, when Γ = [0, 1] n , for standard BEM with a uniform mesh and the same number of degrees of freedom.

Inverse estimates and conditioning
The following inverse estimate follows almost immediately from the wavelet characterisation of the spaces H −t (Γ) for −1 < t < 1 in Theorem 3.1 and Corollary 3.3. In Appendix B we show, by an alternative, lengthier argument, closer to standard arguments based on "bubble functions" (e.g. [20]), that the estimate (86) holds in fact for all t > 0, and for the full range 0 < d < n; moreover, for any T > 0, the constant c t in the estimate can be chosen independently of t for 0 < t ≤ T .
The first claim of the following corollary follows immediately from the above result. Corollary 5.12. Suppose that Γ satisfies the conditions of Theorem 5.6 and that φ and φ N are as defined in Theorem 5.6. Suppose also that −1/2 ≤ s 1 < s 2 < −(n − d)/2 and that φ ∈ H s 2 Γ . Then, for some constant c > 0 independent of h and φ, Recalling from §2.4 that H −s 1 (Γ c ) ⊥ is a unitary realisation of the dual space of H s 1 Γ with respect to the duality pairing ·, · H −s 1 (Γ∞)×H s 1 (Γ∞) , we have that denote the solution of (38) with −g replaced by ϕ. Noting that, by Theorem 2.7, , for some constant C > 0 independent of ϕ. It follows from (80) and (83) that, for some constant c > 0 independent of φ, h, and ζ, Another application of Theorem 5.10, specifically (86), is to prove bounds on the condition number of the matrix in our Galerkin BEM. Let N := #L h and suppose that be the corresponding basis for Y h . Then the Galerkin method using the N -dimensional space V N = Y h leads to the Galerkin matrix A ∈ C N ×N given by (60). The following theorem bounds the 2-norm, · 2 , of this matrix and its inverse. The bound for A 2 is in terms of S 2 , the norm of S as an operator on L 2 (Γ), and recall that t d is defined in (44). The numerical results reported at the end of §6.2 suggest that these bounds are sharp in their dependence on h and d.
Theorem 5.13. Let Γ be an IFS attractor satisfying the OSC with n − 1 < d = dim H (Γ) < n. Then, with V N = Y h as described above, the Galerkin matrix A ∈ C N ×N defined in (60) satisfies If also Γ is disjoint then, for some c > 0, so that also In particular, (90) and (91) hold with c = αc −2 t d , where α is the coercivity constant from (41) and c t d is the constant from (86) when t = t d .

Numerical quadrature and fully discrete error estimates
To evaluate the Galerkin matrix (60) and right-hand side entries (61) we use the quadrature routines from [25]. These are applicable when Γ is a disjoint IFS attractor with n − 1 < d = dim H (Γ) < n, and we assume throughout this section that Γ is of this form.
As in the previous two sections we focus on the case where V N = Y h (with Y h defined as in (74)), for some h ∈ (0, diam(Γ)], and we use the canonical orthonormal basis for V N (i.e., tr * Γ applied to the functions χ m from (75) and (18)) so that A and b are given by (60) and (61). It follows from (60) that the Galerkin matrix entries are the double integrals For i = j the integrand in (93) is continuous, because, by assumption, the components of Γ are disjoint. In that case, to evaluate (93) we use the composite barycentre rule of [25,Defn 3.5].
In our context, this means partitioning the BEM elements Γ m(i) and Γ m(j) , which have diameter approximately equal to h, into a union of self-similar subsets of a possibly smaller approximate diameter 0 < h Q ≤ h, writing the integral in (93) as a sum of integrals over the Cartesian products of these subsets, and approximating these integrals by a one-point barycentre rule. Specifically, we approximate where, for i ∈ 1, . . . , N , describes the partitioning of Γ m(i) and, for n ∈ I N 0 , is the barycentre of Γ n with respect to the measure H d . (Note that x n is not necessarily an element of Γ n .) The similarities can be written as s m (x) = ρ m A m x + v m , for some v m ∈ R n and some orthogonal A m ∈ R n×n (A m = ±1 in the case n = 1), and then, writing n = (n 1 , . . . , n ), we find (see [25,Prop. 3.3]) that the formula that we use for calculation of x n in §6.
The integral of Φ sing over Γ m(i) × Γ m(i) is singular, but, using the self-similarity of Γ and the symmetry and homogeneity properties of Φ sing , namely the fact that, for m = 1, . . . , M , it can be written (see [25,Thm 4.6]) as The integrals appearing in (97) all have continous integrands and can be evaluated using the composite barycentre rule (as in [25,Eqn (47)]), viz.
Let A Q ii denote the resulting approximation of A ii , obtained by combining (95)-(98).
Assuming the data g is sufficiently smooth, the right-hand side entries (61) are regular single integrals and can be evaluated by a single integral version of the composite barycentre rules described for double integrals above, as in [25,Defn 3.1]). Explicitly, we use the approximation Having we solve, instead of (57), the perturbed linear system A Q c Q = b Q . Provided this is uniquely solvable (for which see Corollary 5.17 below), our fully discrete approximation to φ is φ Q N := N j=1 c Q j e j .
Provided ϕ ∈ ( H 1/2 (Γ c )) ⊥ is sufficiently smooth, we approximate J(φ), given by (62), by an approximation J Q (φ Q N ) defined in a similar way to (99). Using (65), we have that The following quadrature error estimates follow from results in [25], specifically [25,Prop. 5 The assumption of hull-disjointness permits analysis of the singular quadrature rules in [25] by Taylor expansion. (The point here is that while the barycentre x m of a component Γ m may not lie in Γ m , it must lie in Hull(Γ m ).) However, numerical experiments suggest that hull-disjointness is not essential for the applicability of the quadrature rules in [25] (see [25, §6]). We also suspect that the estimates (105) and (104) may not be sharp in their h-dependence -see [25,Rem. 5.10] and the discussion around [25, Fig. 7]. In Remark 5.19 we describe modifications to our quadrature rule which may improve its efficiency.
(ii) Suppose that tr Γ ϕ = V | Γ , where V is twice boundedly differentiable in some open neighbourhood of Hull(Γ). (For ϕ given by (63) this holds with V = −v if v is C ∞ in a neighbourhood of Hull(Γ).) Then, for −1/2 ≤ s < −(n − d)/2, there exists a constant C > 0, independent of h, h Q , and V , such that, for all ψ N ∈ V N = Y h , (iii) Suppose that Γ is hull-disjoint. Then there exists a constant C > 0, independent of h and h Q , such that, for i, j = 1, . . . , N , If, further, Γ is homogeneous, then there exists a constant C > 0, independent of h and h Q , such that Proof. In what follows, when applying results from [25] we are taking h, Γ and Γ in [25] to be respectively h Q , Γ m(i) , and Γ m(j) from the current paper. For (i), noting that the first estimate in (102) follows from [25, Thm 3.6(iii)], and the second then follows from the fact that N i=1 µ m(i) = H d (Γ). For (ii), given ψ N ∈ V N let ψ = (( ψ) 1 , . . . , ( ψ) N ) T ∈ C N denote the coefficient vector of its expansion with respect to the basis (e 1 , . . . , e N ). By the orthonormality of the basis (f 1 , . . . , f N ) in L 2 (Γ), the inverse estimate (86) and the isometry property (16), given 0 < t < 1 we can bound Then, to prove (103), using (65), (101), [25, Thm 3.6(iii)], (106), and the Cauchy-Schwarz inequality, For (iii), the first estimates in (104) and (105) follow from [25,Prop 5.2] for the off-diagonal terms and [25,Thms 5.7 & 5.11] for the diagonal terms. The factors of h −(n+1) and h −n come from the fact that, in the notation of [25,Thms 5.7 & 5.11], R Γ,Hull equals diam(Γ) times a constant independent of diam(Γ). We note that, since we are allowing C to be k-dependent, we can use the "k diam(Γ) ≤ c osc " estimates in [25,Thms 5.7 & 5.11], since, given k > 0, the constant c osc can be chosen as large as is required. To derive the second estimate in (104) (a similar argument gives the second estimate in (105)), note that, with · F denoting the Frobenius norm, Thm. 1.14-1.15]; see, e.g., [47] for more recent related results. However, for n > 1 the exact value of the Hausdorff measure of even the simplest IFS attractors is known only for d ≤ 1 (see, e.g., [46]), i.e. for the cases that are not relevant for scattering problems (recall, as discussed above Lemma 4.4, that H −1/2 Γ = {0} and φ = 0 for d ≤ n − 1). Thus, the entries of the Galerkin matrix (60) and right-hand side vector (61) in our Hausdorff BEM can be computed only up to unknown factors that are powers of the unknown measure H d (Γ), so that the same is true for the vector c that is the solution of the Galerkin linear system (57). To be precise, we can't compute c, but we can (even if we don't know H d (Γ)) compute c = (c 1 , . . . , c N ) T , related to c by c = (H d (Γ)) 1/2 c.
Surprisingly, this inability to compute c does not affect in the slightest the applicability of our Hausdorff BEM. Indeed, it is enough to know c (whose value is independent of H d (Γ)) to compute the BEM solution φ N = for every ϕ ∈ H 1/2 (R n ), and the right hand side of the above is independent of the value of H d (Γ).
A simple implementation technique which avoids working with an unknown value for H d (Γ) and produces the correct solution φ N is just to set H d (Γ) = 1 in all calculations. This is equivalent to introducing a "normalised Hausdorff measure" H d (·) := H d (·)/H d (Γ), so that H d (Γ) = 1, and using it throughout in the BEM in place of H d (·).
To study the influence of the quadrature error on the BEM solution we first adapt the first Strang lemma to our setting.
Proposition 5.16. Let Γ be a disjoint IFS attractor with n−1 < d = dim H (Γ) < n. Assume that the unique solution φ of (36) belongs to H s Γ for some −1/2 < s < −(n − d)/2, so that 0 < s + 1/2 < t d . Let V N = Y h so that N = dim(V N ). Let α be the coercivity constant of (41) and c t d the constant in the inverse inequality (86) when t = t d . Assume that A Q ∈ C N ×N and b Q ∈ C N are approximations of the Galerkin matrix A and the right-hand side vector b in (60)-(61) satisfying for some 0 ≤ E A < α/c 2 t d , and E b ≥ 0. Then the perturbed linear system A Q c Q = b Q is invertible and the corresponding solution φ Q N := N j=1 c Q j e j ∈ V N satisfies the error bound where φ N = N j=1 c j e j ∈ V N is the Galerkin solution given by (56), α Q := α − E A c 2 t d , and C > 0 is a constant independent of h, φ, E a , E b , and N .
Proof. As in the proof of Theorem 5.14, given ψ N ∈ V N let ψ ∈ C N denote the coefficient vector of its expansion with respect to the basis (e 1 , . . . , e N ). Let B(ψ N ) := − g, ψ N H 1/2 (Γ∞)×H −1/2 (Γ∞) = ψ H b and denote the perturbed sesquilinear form and antilinear functional by a Q (ξ N , ψ N ) := ψ H A Q ξ and From this and the coercivity of a(·, ·) in (41) follows the coercivity of the perturbed form, and hence the invertibility of A Q : The second bound in (107), combined with (106) for t = t d , gives Strang's first lemma (e.g. and the bound (108) follows on applying Corollary 5.12 with s 2 = s.
By combining Theorem 5.6, Proposition 5.16 and the quadrature error bounds in Theorem 5.14, we can complete the convergence analysis of our fully discrete Hausdorff BEM for a hull-disjoint IFS attractor.
Corollary 5.17. Let Γ be a hull-disjoint IFS attractor with n − 1 < d = dim H (Γ) < n. Assume that the unique solution φ of (36) belongs to H s Γ for some −1/2 < s < −(n − d)/2. Let V N = Y h . Let the entries of the Galerkin matrix A (60) and right-hand side b (61) be approximated with the quadrature formulas outlined in (93)-(99). Let g satisfy the assumptions of Theorem 5.14(i), and suppose that for some sufficiently small C Q > 0, independent of h Q and h. Then the approximated linear system is invertible and the fully discrete solution φ Q N satisfies the error bound for some C > 0 independent of h, φ, and G.
If, further, Γ is homogeneous the above result holds with (111) replaced by the weaker condition Proof. The two conditions (111) and (113) can be written as where p = 0 in the general case, and p = 1 when Γ is also homogeneous. By Theorem 5.14(i) and (iii) our quadrature formulas achieve (107) with for some constant C * independent of both h and h Q . Hence (114) implies that E A is bounded independently of h, and, furthermore, by choosing the constant C Q in (114) to be sufficiently small, one can ensure that E A ≤ α/(2c 2 t d ), so that the discrete system is invertible by Proposition 5.16, with The assumption of (114) also ensures that E b is bounded independently of h, with since h ≤ diam(Γ) and 2d 1+p − s + n−d 2 + 1 > 0 for both p = 0 and p = 1. Then, by (108), for some C > 0 independent of h, φ, and G, and by combining this with (80) we deduce (112).
Under a stronger condition on h Q we can also prove a superconvergence result for the fully discrete approximation J Q (φ Q N ), given by (101), to the linear functional J(φ) given by (62) (cf. [35,Thm. 4.2.18]). The significance of this result for the computation of u(x) and u ∞ (x) is as spelled out in and above Remark 5.7. Corollary 5.18. Suppose that Γ, φ, V N = Y h , and g are defined, and satisfy the same assumptions, as in Corollary 5.17. Let J and J Q be defined by (62) and (101) and suppose that ϕ satisfies the assumptions of Theorem 5.14(ii). Suppose also that ζ ∈ H s Γ , where ζ is the solution of (56) with −g replaced by ϕ, and suppose that for some sufficiently small C Q > 0, independent of h and h Q . Then the approximated linear system is invertible and for some constant C > 0 independent of h, φ, ζ, ϕ, G, and V .
If, further, Γ is homogeneous the above result holds with (117) replaced by the weaker condition Proof. The two conditions (117) and (119) can be written as where p is as in the proof of Corollary 5.17. We first note that, for every C Q > 0, since s > −1/2, (120) implies (114) provided C Q is sufficiently small. Thus, if (120) holds with C Q sufficiently small, then the approximated linear system is invertible by Corollary 5.17. Next we estimate by (81) and Theorem 5.14(ii). We note moreover that φ Q N H s Γ ≤ C 1 ( φ H s Γ + |G| 2,Hull(Γ) ), for some constant C 1 > 0 independent of φ, G, and h, as a consequence of Corollary 5.17 and Proposition 5.11, so that (120) implies that the last term in the above equation is recalling that s > −1/2 so that s( 2 1+p − 1) + 1 + 2 1+p (d + 1/2) + (n − d)/2 > 0 for both p = 0 and p = 1. Further, arguing as in the proof of Corollary 5.17, the conditions (107) are satisfied with E A and E b given by (115), for some constant C * independent of h and h Q . If also (120) holds then 1+p (d+1/2)+1/2+(n−d)/2 |G| 2,Hull(Γ) h s+1/2 , recalling that s < −(n − d)/2, so that (2 − 2 1+p )(−s) + 2 1+p (d + 1/2) + 1/2 + (n − d)/2 > 0 for p = 0, 1. Thus, if C Q is chosen sufficiently small so that E A ≤ α/(2c 2 t d ) for 0 < h ≤ diam(Γ), it follows from Proposition 5.16 that Hull(Γ) )h 2s+1 , for some constant C 2 > 0 independent of h, φ and G. The bound (118) follows by combining the above inequalities. Remark 5.19 (Reduced quadrature). The quadrature rule defined by (94) uses the same maximum mesh width h Q for all off-diagonal elements A ij , i = j. Since the magnitude of the integrand Φ(x, y) and its derivatives blows up as one approaches the diagonal x = y, and decays away from it, it is possible to save computational effort, while maintaining the error bounds (104) and (105), and hence the error estimates (112) and (118). To achieve this we increase the quadrature mesh width (and hence potentially decrease the number of quadrature points) for the computation of A ij when the elements Γ m(i) and Γ m(j) are sufficiently well-separated spatially. In more detail, from [25,Prop 5.2] we have that, for i = j, the error in the quadrature approximation (94), with h Q replaced by a local quadrature mesh width h Q,i,j , satisfies , for some constant C independent of i and j, where Υ(R) := 1 + (kR) n/2+1 R n+1 and R ij := dist(Hull(Γ m(i) ), Hull(Γ m(j) )).
Given h Q , one can therefore maintain the bounds (104) and (105) by replacing h Q in (94) by If the quantities R ij are not known exactly, but satisfy for some known quantities R − ij ≥ 0 and R + ij < ∞, then one can maintain (104) and (105) by replacing h Q in (94) by Here we are using the fact that Υ(R) is a monotonically decreasing function of R, and the max(·, 1) is needed to ensure that h Q,i,j ≥ h Q (so that we are not increasing the computational effort unnecessarily) because the quantity inside the square root may be smaller than 1. (This holds in particular if R − ij = 0, in which case we are interpreting 1/Υ(R − ij ) as 1/Υ(0) = 1/∞ = 0.) In our numerical results in §6 we will compute A Q ij as described above, choosing h Q,i,j according to (123) with and which satisfy (122).

Numerical results
In this section we present numerical results that illustrate our methods and assess the sharpness of our theoretical predictions. The Hausdorff BEM (57) has been implemented in the Julia language [4] and the code is available at https://github.com/AndrewGibbs/IFSintegrals. In all examples described in this section, we consider the scattering of plane waves, i.e. the datum g is as in (31) with u i (x) = e ikϑ·x and |ϑ| = 1. We validate our implementation against a different method in Our implementation uses the quadrature rules described in §5.4. Precisely, we approximate the right-hand side in the linear system and linear functionals of the solution (the scattered field and far field) using the quadrature rules (99) and (100) (113)) and h Q ≤ C Q h 1+d (see (111)) required by our theory, our numerical experiments suggest that, in practice, h Q = C Q h is sufficient to achieve our theoretical convergence rates, even when using the reduced quadrature of Remark 5.19.
We measure the accuracy of our BEM solutions in the H −1/2 Γ norm. These H −1/2 Γ norms are computed by expressing them in terms of a single-layer BIO with wavenumber k = i, which we denote S i , as in [16, Table 1]. Practically, we achieve this by assembling an approximation A i,Q to the Galerkin matrix A i for the operator S i , analogous to (93), with k = i and N ref degrees of freedom, using quadrature approximations analogous to (95)-(98) and (94), choosing h Q the same as for the reference solution, with the reduced quadrature formulae (123)-(125) applied as if k = 1. (While the quadrature convergence analysis in [25, §5] was presented only for k > 0, one can check that the relevant results also extend, mutatis mutandis, to the case k = i.) Next, we view φ Q N as an element of the larger space V N ref and define v as the coefficient vector of Then, arguing as in [16, Table 1] (and see Footnote 16), it follows that We consider first the case where n = 1 (corresponding to scattering in R 2 ) and Γ ⊂ Γ ∞ ∼ = R is a "middle-(1 − 2ρ) Cantor set" screen, for some 0 < ρ < 1/2. Precisely, Γ is the attractor of the disjoint homogeneous IFS with M = 2, a d-set with d = dim H (Γ) = log 2/ log (1/ρ). We denote by φ the Hausdorff-BEM solution "at level ", for ∈ N 0 , i.e. with h = ρ and N = 2 .  In Remark 5.8 we noted that our theoretical analysis suggests we should expect convergence approximately like 2 − /2 for the BEM solution, which is what we observe in the numerical results (for each error convergence plot we also plot C2 − /2 -the black dashed lines -choosing the constant C so that the two curves coincide for the largest value of ). Higher wavenumbers give different pre-asymptotic behaviour (in particular, the approximation is not accurate for small values of ) but do not affect the asymptotic rates for large , see Figure 8. For ρ = 0.5 our Hausdorff BEM coincides with classical piecewise-constant BEM on a uniform mesh applied to the Lipschitz screen that is the unit interval [0, 1]; our implementation uses the special choice (95)-(98) of quadrature rules to evaluate the matrix entries (which is simply the composite midpoint rule for the off-diagonal entries). Running our Hausdorff-BEM code with ρ = 0.5 we observe error curves (not reported here) that are almost identical to those for ρ = 0.49. Hence, at least in this case, the piecewise-constant Hausdorff-BEM approximation of the integral equation on a fractal (the Cantor set with ρ = 0.49) is no less accurate, for the same number of degrees of freedom, than a classical piecewise-constant BEM approximation on an adjacent, more regular set (the interval [0, 1]). We consider next the case where n = 2 (corresponding to scattering in R 3 ) and Γ ⊂ Γ ∞ ∼ = R 2 is a "middle-(1 − 2ρ) Cantor dust" for some 0 < ρ < 1/2, defined by the disjoint homogeneous IFS with M = 4,
In Figure 9 we present results similar to those in Figure 8, albeit for a more restricted range of and a lower value of k in the right-hand panel, due to the increased computational cost associated with the change from M = 2 (Cantor set) to M = 4 (Cantor dust). The Hausdorff-BEM solution φ corresponds to a mesh of N = 4 elements of diameter h = √ 2ρ . The incident wave is the plane wave u i (x) = e ikϑ·x with ϑ = (0, 1, −1)/ √ 2. We use the quadrature rules described in §5.4, with the reduced quadrature of Remark 5.19, and with h Q = ρ 2 h for both k = 0.1 and k = 5. As a result, the off-diagonal entries of the BEM matrix require up to 256 evaluations of Φ.
In Remark 5.8 we noted that our theoretical analysis suggests we should expect convergence approximately like (4ρ) − /2 , whereas the convergence we observe in Figure 9 appears to be faster than this, especially for lower values of ρ. In fact, these numerical results are consistent with the theoretical bounds being sharp, as we now explain. The issue is that the convergence rate (4ρ) − /2 is slow, so the reference solution (with = 8) used in Figure 9 is still quite far from the exact solution φ of the integral equation, which affects the shape of the error plot. Let us assume that the Galerkin error bound (84) holds with = 0 and in fact is exact, i.e.
For ρ = 0.5 our Hausdorff BEM coincides with a piecewise-constant BEM on a uniform mesh applied to the screen [0, 1] 2 ; our implementation uses the quadrature rules (95)-(98) (which reduce to a composite product midpoint rule for the off-diagonal entries). As similarly reported for the Cantor set in §6.1, running our Hausdorff-BEM code with ρ = 0.5 we observe error curves (not reported here) that are almost identical to those for ρ = 0.49.
For this Cantor dust example we also compare numerical results with the theoretical predictions of Theorem 5.13. In Figure 12 we plot the norms of the Hausdorff-BEM Galerkin matrix A and its inverse against (the matrix dimension is N = 4 ) for ρ = 1/3 and different values of k. We observe that A 2 is approximately independent of and h = √ 2ρ , and that, for ≥ 2,

Convergence of the scattered field
We now study the convergence of the Hausdorff-BEM approximation of the field scattered by Γ, for the same scatterers (Cantor sets and dusts), wavenumbers, incident waves and quadrature parameters considered in §6.1-6.2. The near field at x ∈ R n+1 and far-field pattern atx ∈ S n are computed using (100), choosing ϕ = Φ(x, ·) and ϕ = Φ ∞ (x, ·) respectively, using the same h Q values used to construct the associated BEM system.
In Figure 13 we show L ∞ errors in the near-and far-field for scattering by Cantor sets and Cantor dusts, obtained by computing the maximum error over a suitable set of sample points. In more detail, define the parameter N (k) := 10 max(k, 2), which is always an even integer for the values of k in our experiments. For the near-field, when n = 1 we sample at 4N points on the boundary of the square (−1, 2) × (−1.5, 1.5), and when n = 2 we sample at N 2 points on a uniform grid on the square (−1, 2) × (−1, 2) × {−1} (recall that Γ ⊂ [0, 1] × [0, 1] × {0}). For the far-field, when n = 1 we sample at N points on the circle S 1 , and when n = 2 we sample at N 2 × N points on the sphere S 2 , chosen such that the points form a uniform grid in spherical coordinate space [0, π] × [0, 2π]. For Cantor sets, we observe that near-and far-field errors converge to zero with rates 2 − , precisely as predicted by the theory in (84). For Cantor dusts, we observe convergence that is apparently faster than the predicted theory, similar to the observation made in relation to Figure 9, which was explained in §6.2. Applying the same reasoning as in §6.2 to achieve a robust comparison with the theoretical error bounds, in Figure 14 we plot u − u +1 L∞ and u ∞ − u ∞ +1 L∞ against . Both agree with the predicted convergence rate of (ρM ) − .

Comparison against "prefractal BEM"
In this section we compare the approximations produced by our Hausdorff BEM with those produced by the method of [16]. This provides a stronger validation of our Hausdorff BEM than the experiments so far presented, which compare our Hausdorff-BEM approximations against higher-accuracy approximations produced using the same Hausdorff-BEM code. We shall refer to the method of [16] as "prefractal BEM", since it involves applying a standard BEM approximation on a prefractal approximation of Γ, which is the closure of a finite union of disjoint Lipschitz open sets. For brevity we focus only on the case of scattering by Cantor sets. In this case, for the prefractal BEM we take as prefractals the sequence of sets Γ (0) = [0, 1], Γ ( ) := s(Γ ( −1) ), ∈ N, where s is as in (4) for the IFS (126).
For each ∈ N 0 we denote byũ ∞ the prefractal-BEM approximation to the far-field pattern u ∞ , computed on the prefractal Γ ( ) using a standard piecewise-constant Galerkin BEM with one degree of freedom for each connected component of Γ ( ) . Quadrature was carried out using high-order Gauss and product Gauss rules, and quadrature parameters were chosen so that increasing quadrature accuracy did not noticeably change the results presented below. This leads to a total number of degrees of freedom N ( ) = 2 , which is the same number of degrees of freedom used in our Hausdorff BEM at level . As before, we write u ∞ to denote the approximation to the far-field pattern for the Hausdorff BEM, which is computed as described in §6.3, except that now we use smaller values of the quadrature parameter h Q (as detailed below), to ensure that the results presented are not polluted by effects of insufficiently accurate quadrature, so our focus is on the Galerkin error per se.  In Figure 15 we show L ∞ errors in the far-field pattern between the two methods for ρ = 0.1, 1/3 and 0.49, with incident direction ϑ = (1/2, − √ 3/2) and wavenumbers k = 0.1 and k = 50 (with h Q = ρ 6 h for k = 0.1 and h Q = ρ 8 h for k = 50 for the Hausdorff BEM), obtained by computing the maximum error over 300 equally-spaced observation angles in S 1 . For each method we calculate the errors in two different ways, using, firstly, a Hausdorff-BEM reference solution, and, secondly, a prefractal-BEM reference solution; in both cases the reference solution is computed with ref = 13. In all cases, it is clear that both methods are converging to the same solution, validating both methods. Indeed, it appears that both methods converge at the same 2 − rate. However, the results suggest that the relative accuracy of the two methods is dependent on the value of ρ, with the two methods having essentially the same accuracy for ρ = 0.49, the prefractal BEM being more accurate for ρ = 1/3, and the Hausdorff BEM being more accurate for ρ = 0.1. To investigate this further we carried out similar experiments for more values of ρ, plotting the ratio of the errors obtained (measured in the L ∞ and L 2 norms) with = 10 in Figure 16. We observe that: • for ρ between 0.4 and 0.5 the two methods are very similar in accuracy; • for ρ between 0.3 and 0.4 the prefractal BEM appears to be more accurate, significantly so, by a factor > 100, for ρ ≈ 1/3 (this appears to be due entirely to some unexpected enhanced accuracy of the prefractal BEM for ρ ≈ 1/3); • for ρ between 0 and 0.3 the Hausdorff BEM appears to be more accurate, significantly so, by a factor > 1000, for the lowest value of ρ; • the ratio of the errors for the two methods appears to be essentially independent of k for the range of k tested. To illustrate this we have also included in Figure 16 results for k = 30 (computed with h Q = ρ 8 h), alongside the results for k = 0.1 and k = 50; we observe that the results for all three k values are almost identical.
These observations about the accuracy of the two methods (particularly the "spike" in Figure 16 near ρ = 1/3) merit further investigation, but we leave this to future work. Figure 16: Ratio of errors of Hausdorff-BEM (u ∞ 10 ) and prefractal-BEM (ũ ∞ 10 ) approximations of the far-field pattern u ∞ with = 10 for scattering by Cantor sets with k = 0.1, k = 30 and k = 50 and ϑ = (1/2, − √ 3/2), for a range of values of ρ ∈ (0, 1/2), measured in the L ∞ (left panel) and L 2 (right panel) norms on S 1 . A value below 1 indicates that the Hausdorff BEM is more accurate, while a value above 1 indicates that the prefractal BEM is more accurate. In both cases ref = 13.

Non-homogeneous or non-disjoint IFS attractors
We now consider two IFS attractors that are not Cantor sets/dusts, both for n = 2: (i) the "non-homogeneous dust", shown in Figure 17 Both satisfy the OSC (5) and are hence d-sets for the stated values of d. But only (ii) is homogeneous, and only (i) is disjoint and satisfies the assumptions of the Hausdorff-BEM convergence theory of §5.2. Figure 18 shows the H −1/2 Γ norms of the Hausdorff-BEM errors computed against a fixed reference solution (left column) and the norms of the differences between Hausdorff-BEM solutions on consecutive meshes as in Figure 11 Figure 18 appears slightly faster than the theoretical rate h (1−n+d)/2 predicted by Theorem 5.1. As before, this apparent mismatch is due to the limited accuracy of the reference solution; when we plot incremental errors (in the top right panel) we see much clearer agreement with the theory. Example (ii) is not covered by our convergence theory, but if our theory were to extend to this case our predicted convergence rate of h (1−n+d)/2 would evaluate to (2/3) . Similarly to what we observed for (i), the relative errors in the bottom left panel of Figure 18 converge slightly faster than (2/3) , while the increments in the bottom right panel converge approximately in agreement with (2/3) . However, we leave theoretical justification of this empirical observation for future work.

A Besov spaces on d-sets
In this appendix we show that the spaces denoted by B p,q α (Γ), and defined in terms of atoms in [30, section 6], coincide with the spaces denoted by B α p,q,0 (Γ) defined in [8,Def. 6.3], under the assumptions α > 0, 1 ≤ p, q < ∞, and Γ being a d-set, 0 < d < n, preserving Markov's inequality in the sense of [30, §4] (or see [31, p. 34], and note Remark A.6 below). As a consequence (see Corollary A.4) we have that our space H α (Γ), defined in §2.4, coincides with B 2,2 α (Γ) under the same assumptions on α and Γ. 19 We start by rephrasing the above mentioned definition of B p,q α (Γ) in terms closer to the notation used in [7], which we want to use to make the announced connection. In order to do that, we need first to recall the notion of atom as used in [30]. For each ν ∈ N 0 , consider the family of cubes Then, given α, p, d, and K as in Definition A.1 below, an (α, p)-atom a νm associated with Q νm is any a νm ∈ C K (R n ) such that supp a νm ⊂ 3Q νm , |D β a νm (x)| ≤ 2 −ν(α−|β|−d/p) , x ∈ R n , |β| ≤ K, where for a cube Q ⊂ R n and a scalar r > 0 the cube rQ is defined to have the same centre x Q as Q and side length r times that of Q, i.e., rQ := {x Q + r(x − x Q ) : x ∈ Q}.
In [7], atoms are also considered, but the definition is different: for each ν ∈ N 0 , consider the family of cubes Then, given s ∈ R, 0 < p < ∞, K ∈ N 0 , and c ≥ 1, an (s, p) K,0,c -atom a νm associated with Q νm is any continuous function with (classical) partial derivatives up to the order K such that supp a νm ⊂ cQ νm , with convergence in the sense of B

B Inverse estimates
In this appendix we prove that the inverse estimate (86), established for 0 < t < 1 in Theorem 5.10, holds in fact for the extended range 0 < t ≤ J, for any J ∈ N, moreover with a constant in the inverse estimate that is independent of t. Our proof is modelled on standard inverse estimate arguments for negative exponent Sobolev spaces, e.g. [20,Theorem 4.6], with the important difference that it is, of course, impossible to support the smooth "bubble functions" (in the terminology of [20, §4.3]), that are a key tool in the arguments of [20] and in the proof of earlier inverse estimate results, inside an element Γ m ⊂ Γ, as Γ has empty interior. In the case that Γ is a disjoint IFS attractor, it turns out that, to carry through an analogous argument, it is enough to replace the bubble function by a smooth function supported in a carefully chosen neighbourhood of Γ m (the function σ m ∈ C ∞ 0 (R n ) in the proof below). The assumption that Γ is a disjoint IFS attractor implies that Γ is a d-set with 0 < d < n (see Lemma 2.6), but, in contrast to Theorem 5.10, we make no additional constraint on d. We continue to assume in this appendix that h lies in the range (78), i.e. that 0 < h ≤ h 0 := diam(Γ), and to use the notations of the main part of the paper, notably L h , χ m , and Y h defined by (76), (18), and (75), respectively.
Theorem B.1. Suppose that Γ is the attractor of an IFS, satisfying (4), that Γ 1 , ..., Γ M are disjoint (i.e., Γ is disjoint), and that J ∈ N. Then there exists c I > 0 such that, for every ψ h ∈ Y h , Proof. Let ρ min := min and σ m := σ • s −1 m • s −1 m −1 • · · · • s −1 m 1 , and note that, as a consequence of the OSC, the supports of σ m and σ m are disjoint, for m = m . Now suppose that 0 < t ≤ J, and define s in terms of t by (10). For m ∈ L h , noting (15) and (18) Combining this inequality with (138) and (139) and recalling (10) we obtain that and the result follows on recalling that a 2 = ψ h L 2 (Γ) .