Boundary element methods for acoustic scattering by fractal screens

We study boundary element methods for time-harmonic scattering in $\mathbb{R}^n$ ($n=2,3$) by a fractal planar screen, assumed to be a non-empty bounded subset $\Gamma$ of the hyperplane $\Gamma_\infty=\mathbb{R}^{n-1}\times \{0\}$. We consider two distinct cases: (i) $\Gamma$ is a relatively open subset of $\Gamma_\infty$ with fractal boundary (e.g.\ the interior of the Koch snowflake in the case $n=3$); (ii) $\Gamma$ is a compact fractal subset of $\Gamma_\infty$ with empty interior (e.g.\ the Sierpinski triangle in the case $n=3$). In both cases our numerical simulation strategy involves approximating the fractal screen $\Gamma$ by a sequence of smoother"prefractal"screens, for which we compute the scattered field using boundary element methods that discretise the associated first kind boundary integral equations. We prove sufficient conditions on the mesh sizes guaranteeing convergence to the limiting fractal solution, using the framework of Mosco convergence. We also provide numerical examples illustrating our theoretical results.


Introduction
The scattering of acoustic waves by screens (or "cracks" in the elasticity literature) is a classical topic in physics, applied mathematics and scientific computing. The basic scattering problem involves an incident wave propagating in R n (n = 2 or 3), striking a screen Γ , assumed to be a bounded subset (typically, a relatively open subset) of some (n − 1)-dimensional submanifold of R n , and producing a scattered field which radiates outward to infinity. In a homogeneous background medium, the scattering problem can be reformulated as a boundary integral equation (BIE) on the screen, as described in e.g. [35,42,74,75,80], and numerical solutions can then be computed using the boundary element method (BEM), as in e.g. [46,74,75,80]. The classical work cited above has since been extended in many directions, e.g. to the electromagnetic case [13,16], to "multi-screens" (the union of multiple screens intersecting non-trivially) [31], and to hybrid numerical-asymptotic approximation spaces [44] for high-frequency problems.
For simplicity we focus on the case where the screen is flat, and the underlying (n − 1)-dimensional submanifold is the hyperplane Γ ∞ := R n−1 × {0} ⊂ R n . Restricted to this setting, existing studies all assume (either explicitly or implicitly) that the screen Γ ⊂ Γ ∞ is a (relatively) open set with smooth (or piecewise smooth) boundary. In our recent paper [27] we derived well-posed boundary value problem (BVP) and BIE formulations for sound-soft and sound-hard acoustic scattering by arbitrary bounded screens Γ ⊂ Γ ∞ , including cases where Γ or ∂Γ has fractal structure. In this paper we consider the numerical solution of these BVPs/BIEs using the BEM.
Our focus is on scattering by two general classes of fractal screens 1 : (i) Γ is a bounded, relatively open subset of Γ ∞ with fractal boundary, for instance the interior of the Koch snowflake in the case n = 3; (ii) Γ is a compact fractal subset of Γ ∞ with empty relative interior, for instance the Sierpinski triangle in the case n = 3.
In both cases our general approach to analysis and numerical simulation is to approximate the fractal screen Γ by a sequence of smoother "prefractal" screens Γ j , j ∈ N 0 , for which BVP/BIE well-posedness and BEM approximation is classical.
The key question we address in this paper is: Given a fractal screen Γ , how should the prefractals Γ j and the corresponding BEM discretisations be chosen so as to ensure convergence of the numerical solutions on Γ j to the limiting solution on Γ ?
In this paper we focus exclusively on sound-soft screens, on which homogeneous Dirichlet boundary conditions are imposed. Our decision to restrict attention to this case was made partly to make numerical simulations as simple as possible, since one has only to discretise the weakly-singular single-layer boundary integral operator. But the sound-soft case is also particularly interesting from a physical point of view, as there is a strong dependence of the scattering properties on the fractal dimension of the screen (see Proposition 6.2 and the numerical results in § §7. 3-7.4). We leave for future work the application of the techniques developed in this paper to the numerical simulation of the well-posed BVP/BIE formulations for scattering by fractal sound-hard and impedance screens presented in [26,27,43]. (See [26] for simulations of diffraction through fractal apertures in sound-hard screens, equivalent, by Babinet's Principle [14], to the sound-soft screen problem that we focus on in this paper.) Our main results and their novelty. The main focus of this paper is to address the "key question" above, proving results (Theorems 5.2 and 5.3) that specify, for each of the classes (i) and (ii), how to choose a sequence of prefractals and their BEM discretisations so as to achieve convergence of the resulting numerical scheme to the limiting solution on Γ . While BEM simulations have been carried out previously on sequences of prefractals (see Jones et al. [52] in the context of potential theory and Panagiotopulos and Panagouli [63] in elasticity), prior to the results in this paper there does not appear to exist any analysis to justify the convergence of such simulations, that the sequence of numerical solutions converges to the desired limiting solution on Γ .
Our focus throughout the paper is on particular BIEs for sound-soft fractal screens, but we expect the methods and arguments that we introduce to be much more widely applicable. Indeed, our analysis is based on a variational formulation of the BIEs in terms of (complex-valued) continuous sesquilinear forms, which allows the question of prefractal to fractal convergence to be rephrased in terms of the Mosco convergence 2 of the discrete BEM subspaces to the fractional Sobolev space in which the limiting fractal solution lives. The methods that we develop, to reduce proof of convergence of the numerical solution to the Mosco convergence of the BEM subspaces (Lemma 2.5), and to prove Mosco convergence of the BEM subspaces, are potentially widely applicable to Galerkin discretisations of other integral or differential equations posed on rough (not necessarily fractal) domains that are approximated by more regular sets. For example, the proof of Theorem 5.3 depends only on a characterisation of Mosco convergence (Lemma 2.4), quantitative bounds on the norms of mollification operators on scales of Sobolev spaces (Appendix B), and a quantitative extension of standard piecewise constant finite element approximation theory (Lemma A.1), all of which should be widely applicable.
A feature of our numerical analysis based on Mosco convergence is that our discrete BEM subspaces need not be subspaces of the Hilbert space in which the solution on Γ lies; indeed, this is crucial whenever the prefractals are not subsets of the limiting fractal Γ . Our Lemma 2.5, which applies in such cases (and is proved in slightly more generality than we need for the results in §5, anticipating wider application) can be seen as a replacement in this circumstance for the standard Céa lemma and its generalisation to compact perturbations (see, e.g. [73,Thms. 8.10,8.11]). 2 Mosco convergence (Definition 2.2 below) is a standard notion in the study of variational inequalities, closely related to Gamma convergence. Introduced by U. Mosco in [58] (almost exactly 50 years ago), it has been applied by a number of authors to the study of PDEs on sequences of domains, see e.g. [5,15,34,57] and the references therein. It has been mainly used (e.g. the references just cited) in relation to convergence at a continuous rather than a discrete level (i.e. in the context of mathematical analysis rather than numerical analysis). But it is also relevant to proving convergence of numerical methods, as was illustrated already in [59,Chapt. 3] in the context of numerical methods for variational inequalities.
While our main results, Theorems 5.2 and 5.3, are quite general in terms of geometry of the screen Γ and its approximating prefractal sequence, we pay particular attention in our theory, examples, and numerical simulations to cases where Γ (or the boundary of Γ ) is the fixed point of an iterated function system (IFS) (e.g. Corollaries 4.6 and 5.8). In particular our examples in §6 and §7 include the cases where Γ is (when n = 2) a Cantor set and (when n = 3) a Cantor dust, the Sierpinski triangle, or the interior of a Koch snowflake. These are cases where the BIEs we wish to solve are posed either on fractal sets or on rough domains with fractal boundaries. Unsurprisingly, subtle and interesting properties of fractional Sobolev spaces and integral operators on rough sets, explored recently in [17,25,28,29,45], are crucial to our arguments throughout.
A novel feature of our BEM and its analysis is that convergence can be achieved in regimes where each boundary element contains many disjoint components of a prefractal (e.g. Corollary 5.8(ii) and Figure 1). To justify this we need an extension, that applies in such cases, of standard piecewise constant approximation theory in scales of Sobolev spaces with explicit constants. We supply this in Lemma A.1.
Applications and motivations. Wave scattering by fractal structures is relevant for numerous applications, since fractals provide a natural mathematical framework for modelling the multiscale roughness of many natural and man-made scatterers. We highlight in particular the propagation of acoustic and electromagnetic waves in dendritic structures like the human lung in medical science [1,51], and the scattering of electromagnetic waves by snowflakes, ice crystals and other atmospheric particles in climate science [9,22,72,78]. But particularly close to the fractal screen scattering problems that we focus on in this paper are configurations arising in the design of fractal antennas for electromagnetic wave transmission/reception (see e.g. [39,65,81]) and fractal piezoelectric ultrasound transducers (see e.g. [2,3,38,61]) (fractal structures being attractive in these contexts because of the possibility of wideband performance), and configurations that arise in the study of fractal aperture problems in laser physics [30,47,50,79]. The current study into acoustic scattering by fractal screens represents a first step towards the rigorous numerical analysis of integral equation methods for the study of such challenging problems involving fractal scatterers.
Related literature. One of the first studies of wave scattering by fractals appears to be M. Berry's 1979 paper [12] on scattering by "random phase screens", in which Berry coins the term "diffractal" to describe waves that have undergone interactions with fractal structures. The difficult problem of studying high frequency asymptotics for fractal scattering problems was investigated by Sleeman and Hua in [48,70]. Concerning the study of PDE problems on fractal domains more generally, U. Mosco notes in [60] that ". . . introducing fractal constructions into the classic theory of PDEs opens a vast new field of study, both theoretically and numerically", but also that "this new field has been only scratched". In addition to papers cited above, research that has started to explore this field of study includes work on finite element method approximations of heat transmission across fractal interfaces [6,18,20,21,55] and H 1 extension problems [36]; Dirichlet-to-Neumann (Poincaré-Steklov) operators on domains with fractal boundaries [4]; and finite difference [62] and conformal mapping [7] approaches to the computation of Laplace eigenfunctions on fractal domains.
Structure of the paper. In §2 we collect some important Hilbert and Sobolev space results that will be used throughout the paper. The main new result here is Lemma 2.5, which proves convergence of solutions of variational problems for "compactly perturbed coercive" sesquilinear forms on a Mosco-convergent sequence of closed subspaces. The compactly perturbed coercive setting is more general than we need for the particular problem under consideration in this paper, since for flat soundsoft screens the first-kind formulation is coercive (strongly elliptic) [25]. However, it is included here to lay the foundations for future investigations into other, closely related problems, such as curved sound-soft screens (as studied e.g. in [75,80]), and impedance screens (as studied e.g. in [8,11,43,54]). In §3 we state well-posed BVPs and BIEs for open and closed screens, refining results from [27], and in §4 we prove convergence of solutions on prefractal screens to solutions on limiting fractal screens using the Mosco framework; in particular we prove for the first time convergence in cases where the prefractal sequence is not monotonic, and it holds neither that Γ ⊂ Γ j for all j, nor that Γ j ⊂ Γ for all j. In §5 we study numerical discretisations based on piecewise constant BEM approximations on prefractals, determining conditions under which the BEM solution converges to the limiting fractal solution in the joint limit of prefractal and mesh refinement. In §6 we present examples of the kind of fractal screens we have in mind, and in §7 we provide numerical results which illustrate our theoretical predictions through a number of concrete examples. We also include in this section preliminary numerical investigations into physical questions such as how the fractal dimension of a screen affects the magnitude of the resulting scattered field. In §8 we make a brief conclusion and list some of the many intriguing open problems.

Dual space realisations
We say that a linear isomorphism between Hilbert spaces is unitary if it preserves the inner product (equivalently, if it is isometric [33,Prop. 5.2]). If H is a complex Hilbert space (as all the Hilbert spaces in this paper are) by its dual space H * we mean, following Kato [53], the space of anti-linear continuous functionals on H. It is often convenient to identify H * , itself a Hilbert space with the standard induced operator norm, with some other Hilbert space H, termed a realisation of H * . If ·, · is a continuous sesquilinear form on H × H, and the mapping taking φ ∈ H to φ ∈ H * , given by φ (ψ) = φ, ψ , ψ ∈ H, is a unitary isomorphism, then we say that H is a unitary realisation of H * with associated duality pairing ·, · . If W ⊂ H is a closed subspace, then a unitary realisation of W * is provided by the following simple but important result, e.g. [ We say that a and A are compactly perturbed coercive if a = a 0 + a 1 with a 0 coercive and a 1 compact; equivalently, A = A 0 + A 1 with A 0 coercive and A 1 compact. 3 Let W ⊂ H be a closed subspace of H and let W := (W a,H ) ⊥ ⊂ H be the unitary realisation of W * provided by Lemma 2.1. We say that a is invertible on W if, for every f ∈ W, the variational problem: find u W ∈ W such that has exactly one solution u W ∈ W . This holds (e.g. [49,Thm. 2.15]) if and only if In terms of the associated operator A, (2) can be written equivalently as where P W is orthogonal projection in H onto W, so that a is invertible on W if and only if P W A| W is invertible, in which case (P W A| W ) −1 = β −1 . If a is coercive then a is invertible on W by the Lax-Milgram lemma and β ≥ α. More generally, if a is compactly perturbed coercive, then a is invertible on W if and only if it is injective, meaning that the problem (2) has at most one solution u W ∈ W for every f ∈ W.

Mosco convergence
We now consider the problem of approximating the solution of the variational problem (2) by the solutions of variational problems posed on a sequence of closed subspaces (W j ) ∞ j=1 ⊂ H. We say that a is uniformly invertible on such a sequence (W j ) ∞ j=1 if a is invertible on W j for all j ∈ N and the inverses are uniformly bounded, meaning that, for some constant C > 0 and all j ∈ N and f j ∈ W j := (W a,H j ) ⊥ , where u W j is the unique solution of (2) with W replaced by W j and f by f j . Equivalently, a is uniformly invertible on (W j ) if P W j A| W j is invertible for j ∈ N and Roughly speaking, given a variational problem (2), to ensure that the corresponding solutions on (W j ) ∞ j=1 converge to the solution in W we require that a is sufficiently "well-behaved" and that W j approximates W increasingly well as j → ∞ in an appropriate sense. The precise requirement on W j is that it converges to W in the Mosco sense (Lemma 2.7 below), and Lemma 2.5 makes clear that a being compactly perturbed coercive is sufficiently "well-behaved". The following definition of Mosco convergence is precisely the notion of set convergence for convex sets introduced in [58, Definition 1.1], except that we specialise here from the general Banach space setting to the Hilbert space case, and our convex sets are specifically closed linear subspaces (as, for example, in [57] (i) For every w ∈ W and j ∈ N there exists w j ∈ W j such that w j − w H → 0.
(ii) If (W j m ) is a subsequence of (W j ), w m ∈ W j m , for m ∈ N, and w m w as m → ∞, then w ∈ W .
The following "sandwich lemma" is a trivial consequence of Definition 2.2.
j for each j ∈ N, and both W + j and W − j Mosco-converge to some closed subspace W of H, then W j also Mosco-converges to W .
The following lemma will also be useful. Again the proof is straightforward. Lemma 2.4. Let W j and W be closed subspaces of H. To prove that W j M −→ W it suffices to show that (i) there exists a dense subspace W ⊂ W such that for every w ∈ W and j ∈ N there exists w j ∈ W j such that w j − w H → 0, and (ii) there exists a sequence of closed subspaces W + j of H such that W j ⊂ W + j for all j ∈ N and W + j M −→ W .
Mosco convergence was introduced to study convergence of approximate solutions to variational inequalities. The following lemma, which applies to variational equalities, appears to be new, but its proof, if specialised to the coercive case, has something of the flavour of the original arguments of Mosco [58,Thm. A and its Corollary]. Indeed, in the case that H is a real Hilbert space and a is coercive, this lemma is a corollary of the results in [58], since variational inequalities on linear subspaces of real Hilbert spaces are in fact equalities.
Lemma 2.5. Let W ⊂ H and W j ⊂ H, for j ∈ N, be closed subspaces such that W j Mosco-converges to W as j → ∞. Let a be compactly perturbed coercive and invertible on W . Then there exists J ∈ N such that a is uniformly invertible on W j for j ≥ J. Further, if, for some f ∈ H, u W denotes the solution to (2) and, for j ≥ J, u W j denotes the solution to (2) with W replaced by W j , then Proof. We show first that, for some J ∈ N, a is uniformly invertible on W j for j ≥ J. Suppose first that a is not invertible on W j for all sufficiently large j, in which case neither is it injective on W j . Then there exists a subsequence of (W j ), which we will denote again by (W j ), and v j ∈ W j with v j H = 1, such that If on the other hand a is invertible on W j for all sufficiently large j but is not uniformly invertible then there exists a subsequence of (W j ), which we will denote again by (W j ), and v j ∈ W j with v j H = 1 such that In both of these cases as v j H = 1 is bounded we can extract a subsequence that is weakly convergent to some v ∈ H. Denoting the subsequence again by (v j ), we have that v j ∈ W j and v ∈ W by (ii) in Definition 2.2. Further, by (i) in Definition 2.2, for all w ∈ W , there exists a sequence (w j ) ⊂ H with w j ∈ W j such that w j − w H → 0. Thus, and by (6) or (7), But also a(v j , w) → a(v, w). Thus a(v, w) = 0 for all w ∈ W , so that v = 0 as a is invertible on W . So v j v = 0 and, by (6) or (7), a(v j , v j ) → 0. Further, recalling that a = a 0 + a 1 with a 0 coercive and a 1 compact, we have also a 1 (v j , v j ) → 0 as a 1 is compact. But this implies that a 0 (v j , v j ) → 0, which contradicts that a 0 is coercive. Thus, for some J ∈ N, a is invertible on W j for j ≥ J and is uniformly bounded.
Thus, for j ≥ J, u W j is well-defined and (u W j ) ∞ j=J is bounded and so has a weakly convergent subsequence, converging to a limit u * , and u * ∈ W by (ii) in Definition 2.2. Further, by (i) in Definition 2.2, for all w ∈ W there exists w j ∈ W j such that w j − w H → 0, and (2) gives as j → ∞ through that subsequence. Thus a(u W , w) = a(u * , w), for all w ∈ W , so that u * = u W by the invertibility of a on W . By the same argument every subsequence of (u W j ) ∞ j=J has a subsequence converging weakly to u W , so that (u W j ) ∞ j=J converges weakly to u W . Finally, we see that as j → ∞, by the weak convergence of (u W j ) ∞ j=J and (2). Since a 1 is compact, Remark 2.6. The statement of Lemma 2.5 can be strengthened if additional assumptions are made on a, (W j ) ∞ j=1 and W .
(i) If a is coercive then one can take J = 1 (since a coercive sesquilinear form is automatically invertible on every subspace). In the special case when (5) holds, this was noted already in [29,Lem. 2.4]. (ii) If W j ⊂ W for each j ∈ N then (ii) in Definition 2.2 holds automatically and quasi-optimality holds asymptotically [73,, meaning that, for some M > 0 and J ∈ N, a is invertible on W j for j ≥ J, and Furthermore, if a is coercive then, by Céa's lemma, (8) holds with J = 1 and M = C/α, where C and α are the continuity and coercivity constants for a.

Sobolev spaces and trace operators
Our notation follows that of [56] and [29]. Let m ∈ N. For a subset E ⊂ R m we denote its complement E c := R m \ E, its closure E and its interior E • . We denote the Hausdorff (fractal) dimension of E by dim H E (see, e.g., [76, §I.2]). We say that a non-empty closed set F ⊂ R m is a d-set for some 0 ≤ d ≤ m, if it is "uniformly d-dimensional", more precisely if there exist c 1 , c 2 > 0 such that . For m = 1 these definitions coincide: we interpret them both to mean that Ω is a countable union of open intervals whose closures are disjoint and whose endpoints have no limit points. For s ∈ R, let H s (R m ) denote the Hilbert space of tempered distributions whose Fourier transforms (defined for ξ ∈ R m asû(ξ) := In particular, H 0 (R m ) = L 2 (R m ) with equal norms. For the dual space of H s (R m ) we have the unitary realisation (H s (R m )) * ∼ = H −s (R m ), with duality pairing which coincides with the L 2 (R m ) inner product when both u and v are in L 2 (R m ). Given a closed set F ⊂ R m , we define and given a non-empty open set Ω ⊂ R m , we define Although H s (Ω) is a space of distributions on Ω, it can be naturally identified with a space of distributions on R m , namely (H s Ω c ) ⊥ ⊂ H s (R m ), where ⊥ denotes orthogonal complement in H s (R m ), with the restriction operator | Ω : (H s Ω c ) ⊥ → H s (Ω) providing a unitary isomorphism between the two spaces.
Regarding duality, for arbitrary F ⊂ R m closed and Ω ⊂ R m open, we can unitarily realise the dual spaces of H s F and H s (Ω) by certain closed subspaces of H −s (R m ), with the duality pairing inherited from (9). Precisely, by Lemma 2.1,  (9). An alternative, and more widely-known unitary realisation of ( H s (Ω)) * (also valid for arbitrary open Ω ⊂ R m , see [29,Thm. 3.3]) is where U ∈ H −s (R m ) is any extension of u ∈ H −s (Ω) with U | Ω = u. That (11) and (12) are both unitary realisations of ( H s (Ω)) * is consistent with the fact that is a unitary isomorphism, as mentioned above. In the context of our screen scattering problem, we define Sobolev spaces on the hyperplane Γ ∞ = R n−1 × {0} by associating Γ ∞ with R n−1 and setting H s (Γ ∞ ) := H s (R n−1 ), for s ∈ R, which we shall frequently abbreviate to simply In the exterior domain D := R n \ Γ we work with Sobolev spaces defined via weak derivatives. Given a non-empty open Ω ⊂ R n , let W 1 (Ω) := {u ∈ L 2 (Ω) : ∇u ∈ L 2 (Ω)} and let W 1,loc (Ω) denote the "local" space in which square integrability of u and ∇u is required only on bounded subsets of Ω. We note that, typically, H 1 (D) W 1 (D), since the restriction space inherits from H 1 (R n ) a requirement of (weak) continuity across Γ . We define U + := {(x 1 , . . . , x n ) ∈ R n : x n > 0} and U − := R n \ U + , and adopt the convention that the unit normal vector n on Γ ∞ points into U + . From the half spaces U ± to the hyperplane Γ ∞ we define the standard trace operators γ ± : . We shall frequently abuse notation and apply γ ± and ∂ ± n to elements u of the local space W 1,loc (D), assuming implicitly that u has been pre-multiplied by a cutoff φ ∈ C ∞ 0 (R n ) satisfying φ = 1 in some neighbourhood of Γ , and restricted to U ± as appropriate; for example, γ + u should be interpreted as γ + ((uφ)| U + ).

Boundary value problems and boundary integral equations
Given a screen Γ ⊂ Γ ∞ (a bounded subset of Γ ∞ ), and an incident wave u i ∈ H 1,loc (R n ) (for instance, a plane wave u i (x) := e ikd·x with d a unit direction vector), we seek a scattered acoustic field u satisfying the Helmholtz equation in D = R n \ Γ , the Sommerfeld radiation condition and the Dirichlet boundary condition To formulate a well-posed BVP, one needs to be more precise about the sense in which the boundary condition (15) holds. A detailed investigation into this issue for general bounded subsets Γ ⊂ Γ ∞ was carried out in [27]. Here we apply the results of [27] to describe well-posed BVP formulations for the two types of screen mentioned in the Introduction, namely (i) bounded, relatively open subsets of Γ ∞ , and (ii) compact subsets of Γ ∞ , possibly with empty relative interior. From now on, for brevity we shall omit the words "relative" and "relatively" when discussing relatively open subsets of Γ ∞ and relative complements, boundaries and interiors of subsets of Γ ∞ .

Well-posed BVPs and BIEs for bounded open screens
Let Γ be a non-empty bounded open subset of Γ ∞ . Then we can formulate the scattering BVP by imposing the Dirichlet boundary condition (15) by restriction to Γ (denoting this problem as D(Γ ) r , D for Dirichlet, r for restriction).  (14), and A well-posedness result for this formulation is provided in Theorem 3.2, which is proved in [27, Thm. 6.2(a)]. Before stating it we need some more notation. For Γ ⊂ Γ ∞ non-empty, bounded and open let S Γ : 0 (k|x − y|) (n = 2), and S r Γ : H −1/2 (Γ ) → H 1/2 (Γ ) the single layer boundary integral operator (BIO), the bounded linear operator defined by The key condition for the well-posedness in Theorem 3.2 is H −1/2 (Γ ) = H −1/2 Γ ; the next proposition gives sufficient conditions on Γ for this to hold.

Proposition 3.4. Each of the following are sufficient for
(ii) Γ is C 0 except at a set of countably many points P ⊂ ∂Γ such that P has only finitely many limit points [29,Thm. 3.24]; In §6 we shall combine Theorem 3.2 and Proposition 3.4 to obtain well-posedness results for three-dimensional scattering by certain generalisations of the classical Koch snowflake, as immediate corollaries of the recently established thickness results in [17] (which build on earlier results in [77, Prop. 3.8(iii)]). We shall also deduce well-posedness results for scattering by the standard prefractal approximations to various well-known fractals including the Koch snowflake (and its generalisations), the Sierpinski triangle, and the Cantor dust. In all these cases the standard prefractals are either C 0 or C 0 except at a finite set of points.
Recalling from §2.4 that | Γ : (H is a unitary isomorphism, we note that the problem D(Γ ) r can be equivalently stated with the boundary condition (15) imposed by orthogonal projection. (See Remark 3.6 below for an explanation of why this makes sense physically.) It is instructive to write down this equivalent formulation explicitly, since we will adopt a similar viewpoint when defining BVPs for scattering by compact screens in §3.2. In the following let P Γ : Γ c ) ⊥ (specifically g := −P Γ γ ± u i for the scattering problem), find u ∈ C 2 (D) ∩ W 1,loc (D) satisfying (13) in D, (14), and the boundary condition Remark 3.6. To understand why the boundary condition (20) makes sense in the scattering problem, let u t := u + u i be the total field (the sum of the scattered and incident fields), and consider the traces γ ± u t of u t on Γ ∞ ⊃ Γ . According to formulation D(Γ ) r these traces vanish on Γ , so their supports lie in the complement Since the restriction operator | Γ : (H of both problems can be represented as (18) is the common unique solution of the BIE (19) and the BIE 3.2 Well-posed BVPs and BIEs for compact screens Now let Γ ⊂ Γ ∞ be compact. In particular we have in mind the case where Γ has empty interior, in which case it is not possible to impose the boundary condition (15) by restriction. But, inspired by Definition 3.5 and Proposition 3.7, we can impose (15) by an appropriate orthogonal projection; we justify this in Remark 3.9 below. Extending our existing notation, for compact Γ let P Γ denote the orthogonal projection P Γ : Definition 3.8 (Problem D(Γ ) for compact screens). Let Γ ⊂ Γ ∞ be non-empty and compact. Given (14), and Remark 3.9. We can justify the formulation in Definition 3.8, in particular (22), by relating it to a more familiar formulation of the scattering problem for compact screens [27, Def. 3.1] that replaces (22) with the requirement that . This formulation is well-posed [27, Thm 3.1] when Γ is any compact subset of R n . In the particular case when Γ ⊂ Γ ∞ is a screen it is easy to see that u t ∈ W 1,loc Before stating a well-posedness result for this formulation (which we do in Theorem 3.10), we need some more notation. Given Γ ⊂ Γ ∞ compact and > 0, Remark 3.11. Suppose that Γ ⊂ Γ ∞ is non-empty, bounded and open, in which case Γ is compact, and suppose also that Then we have a choice of well-posed formulations for the scattering problem, potentially with different solutions: problem D(Γ ) (see Definition 3.5) with g := −P Γ γ ± u i (equivalently, [29,Lem. 3.26] and the proof of [17,Lem. 4.15]), and from this it follows that S Γ = S Γ , P Γ = P Γ , and S Γ = S Γ = | −1 Γ • S r Γ . So the two problems D(Γ ) and D(Γ ) are equivalent, sharing the same unique solution.
It is natural to ask whether, in the case where Γ ⊂ Γ ∞ is compact with empty interior, the screen scatters waves at all. This question was answered in [27]. Let Γ ⊂ Γ ∞ be non-empty and compact.
The question of whether H s K = {0} for given s ∈ R and compact K ⊂ R m was investigated in detail in [45]. For sets of Lebesgue measure zero, the Hausdorff [17]. The following lemma collects the results from [17,45]

Variational formulations of the BIEs
To state and analyse Galerkin methods for the BIE formulations, in particular to use the Mosco convergence theory of §2.3, we need variational formulations.
When Γ is a bounded open set and H −1/2 (Γ ) = H −1/2 Γ , or when Γ is compact, we have written down, in Proposition 3.7 and Theorem 3.10, BIEs that are equivalent to the corresponding BVP formulation D(Γ ). In each case these take the form (21) in the respective cases. In each case S Γ : V (Γ ) → V (Γ ) * is a bounded linear operator, a version of the single-layer potential operator.
As noted in §2.2, S Γ has an associated sesquilinear form a Γ (·, ·) defined by The second equality in (24) holds since S Γ = P Γ γ ± S Γ , where P Γ is orthogonal projection onto V (Γ ) * , and since (V (Γ ) * ) ⊥ is the annihilator of V (Γ ), as noted below (11). A consequence of (24) and the definition of i.e. a Γ (·, ·) is the restriction to V (Γ ) of a Γ † (·, ·). Since we can choose Γ † to be as smooth as we wish, e.g. C ∞ , or indeed just an open disk, we see that, even when Γ is fractal or has fractal boundary, the sesquilinear forms we have to deal with are no more complicated than in the case when Γ is a disk. A further consequence of (24) is that, in the case when Γ is bounded and open, (26) so that, where S r Γ is the more familiar screen single-layer potential operator defined by (17), a Γ (·, ·) is also the sesquilinear form associated to S r Γ . Whether Γ is bounded and open or is compact, the sesquilinear form a Γ (·, ·) is continuous. Further, as a consequence of our assumption that the screen is flat For the bounded open case this is hinted at in [42,Rem. 6] and proved rigorously in [25], the latter reference also detailing the wavenumber-dependence of the continuity and coercivity constants. That coercivity of a Γ (·, ·) holds also for every compact Γ is a simple consequence of (25), since coercivity implies coercivity on every closed subspace.
As noted in the general Hilbert space setting in §2.2 (see (2) and (3)), the variational problem: is equivalent to the BIE S Γ φ = −g. The duality pairing on the right hand side can be written equivalently as Since a Γ (·, ·) is coercive, (27) is well-posed by the Lax-Milgram lemma.

Prefractal to fractal convergence
Now suppose we want to study a sequence of problems on a sequence of screens (Γ j ) j∈N 0 (each non-empty and either bounded and open or compact) approximating a limiting screen Γ (again, non-empty and either bounded and open or compact). Assuming that the Γ j are uniformly bounded, let Then by the continuity and coercivity of the are well-posed. The following theorem follows from Lemma 2.5, combined with the continuity of S Γ † : (16)). (28) and (29) respectively Definition 4.2 (BVP convergence). Suppose that Γ is non-empty and either bounded and open or compact, and that the sequence (Γ j ) j∈N 0 is uniformly bounded, and that each Γ j is non-empty and either bounded and open or compact. Then we (23).
The rationale behind this terminology is that, with these conditions on Γ and and φ j are the solutions to (28) and (29), respectively, and that S Γ j φ j → S Γ φ in W 1,loc (R n ), where (noting (25) and the equivalence of (26) with the BIE) S Γ φ and S Γ j φ j are the unique solutions to the BVPs D(Γ ) and D(Γ j ), with data g := P Γ g † and g j := P Γ j g † . In particular, S Γ φ and S Γ j φ j are the scattered fields for scattering of the incident field u i by Γ and Γ j , respectively, provided g † := −P Γ † γ ± u i . Sufficient conditions guaranteeing BVP convergence are given in the following proposition, which follows trivially from (4) where Γ and Γ j are non-empty, bounded and open with Γ j ⊂ Γ j+1 , j ∈ N 0 , and Γ = j∈N 0 Γ j ; where Γ and Γ j are non-empty, bounded and open with H −1/2 (Γ ) = H −1/2 Γ , and there exist Γ − j non-empty, bounded and open and Γ + j non-empty and compact such that: , that applies in more subtle cases where neither Γ j ⊂ Γ nor Γ ⊂ Γ j , is new. We present an example of this type (the "square snowflake") in §6.5 below.

Fractals that are attractors of iterated function systems
While Proposition 4.3 (ii), and its BEM version in Proposition 5.3 below, apply to more general compact Γ , our main motivation for these results is the case when Γ is fractal. An important fractal class (e.g. [37,Chap. 9]) is the set of fractals obtained as attractors of an iterated function system (IFS) {s 1 , s 2 , . . . , s ν }. Here ν ≥ 2 and each s m : R n−1 → R n−1 is a contraction, meaning that for some c ∈ (0, 1). The attractor of the IFS is the unique non-empty compact set Γ satisfying That (31) has a unique fixed point follows from the contraction mapping theorem since s is a contraction on the set of compact subsets of R n−1 , a complete metric space equipped with the standard Hausdorff metric, e.g. [37, Thm. 9.1 and its proof]. If Γ 0 is any non-empty compact set then the sequence Γ j defined by converges in the Hausdorff metric to Γ . In particular, if Γ 0 is such that s(Γ 0 ) ⊂ Γ 0 then [37, Thm. 9.1] In the case that Γ is a fractal or where Γ is not fractal but has a fractal boundary it is common to refer to Γ j as a sequence of prefractals.
The following is an obvious corollary of the above observations, Theorem 4.1, and Proposition 4.3(ii).
Remark 4.7. With Γ j defined by (32), it holds for any non-empty compact Γ 0 that Γ j → Γ in the Hausdorff metric. However, if s(Γ 0 ) ⊂ Γ 0 it may or may not hold that φ j → φ as j → ∞. In particular, if: i) Γ 0 is a countable set; or ii) n = 3 and dim H Γ 0 < 1; then Γ j defined by (32) is also countable or has dim H Γ j < 1, respectively. (The latter case is a consequence of [37,Prop. 2.3].) In such cases it follows from Lemma 3.13 that H

Boundary element methods and their convergence
In this section we propose Galerkin boundary element methods, based on piecewise constant approximations, for solving the screen scattering problems of §3 and prove their convergence. These methods are based on discretisation of a sequence of more regular screens (Γ j ) j∈N 0 ; as in the previous section these converge in an appropriate sense to a limiting screen Γ ⊂ Γ ∞ for which we wish to compute the solution to the scattering problem D(Γ ) of Definition 3.5 or 3.8.
In particular we have in mind cases where Γ is a compact set with fractal dimension < n − 1, and cases where Γ is a bounded open set with fractal boundary. We prove convergence results that apply in both of these cases, and our examples in §6 and our numerical results in §7 are of these types. In each of our convergence results, Γ j is a sequence of open sets, divided into a mesh of elements. Appropriately, given that we are approximating on very rough domains, the constraints on the elements are mild compared to conventional BEM results. In particular our elements need not be convex or even connected.
In more detail, we assume each Γ j is a non-empty bounded open set. On each Γ j we construct a pre-convex mesh M j = {T j,1 , T j,2 , . . . , T j,N j } in the sense of Appendix A, meaning that T j,l ⊂ Γ j is non-empty and open for l = 1, . . . , N j , the convex hulls of T j,l and T j,l are disjoint for l = l, ∂T j,l has zero (n − 1)dimensional Lebesgue measure, and Γ j is the interior of the union of the closures of the T j,l , i.e.
We call h j := max l∈{1,...,N j } diam(T j,l ) the mesh size and T j,1 , T j,2 , . . . , T j,N j the elements of the mesh.
Our Galerkin boundary element method (BEM) is to solve the variational problem (29) with V j chosen to be the N j -dimensional space of piecewise constant functions on the mesh M j , which we denote by V h j . It follows from (25), (26), and the comment following (9), that the BEM solution φ h j is defined explicitly by where (·, ·) L 2 is the inner product on L 2 = L 2 (Γ ∞ ) = L 2 (R n−1 ). Moreover, when we are solving the scattering problem, g † = −P Γ † γ ± u i , and (34) can be written as Definition 5.1 (BEM convergence). Let the discrete approximation space V h j be defined as above, and let the Sobolev space V (Γ ) be defined as in (23).
Mosco-converges to V (Γ ) then we say that "BEM convergence holds". In this case, it follows by Theorem 4.1 that φ h j → φ in H −1/2 , where φ h j and φ are the solutions to (34) and (28), respectively, and

Compact screens
When the limiting screen Γ is compact and H −1/2 (Γ • ) = H −1/2 Γ , neither Theorem 5.2 nor its method of analysis can be applied. In particular, if H −1/2 Γ = {0} and Γ has empty interior (Γ • = ∅) then it is clearly impossible to approximate a limiting non-trivial integral equation solution v ∈ H −1/2 Γ by a sequence of elements of C ∞ 0 (R n−1 ) supported inside Γ , since no such non-trivial functions exist. In the following theorem we address this case. The proof relies on mollification arguments (see Appendix B) to obtain smooth approximations to v to which we can apply the BEM approximation theory (Lemma A.1). This produces approximating smooth functions whose support is strictly larger than that of v. This introduces a constraint on the sequence Γ j to which the analysis applies. In particular, each Γ j must contain Γ ( ) := {x ∈ Γ ∞ : dist(x, Γ ) < }, the -neighbourhood of Γ , for some carefully chosen j-dependent > 0. As is the case throughout this section, our approximation space on Γ j remains the space V h j of piecewise constants on a mesh M j .
, and (Note that we are including t = −1/2 as a possibility here, in which case density holds trivially.) For each j ∈ N defineṽ j := ψ j /2 * v to be the mollification defined in Appendix B. Thenṽ j ∈ C ∞ 0 (Γ j ) (sinceṽ j is smooth and suppṽ j ⊂ Γ ( j /2) ⊂ Γ ( j ) ⊂ Γ j ) and ṽ j − v H −1/2 (Γ ∞ ) → 0 as j → ∞, since j → 0. It remains to show that there exists v j ∈ V h j such that v j −ṽ j H −1/2 (Γ ∞ ) → 0 as j → ∞. For this we define v j to be the orthogonal projection in Hence pairwise disjoint, congruent subsets of Γ ∞ that tile Γ ∞ in the sense that Γ ∞ is the closure of ∞ n=1 S j,n . Let h j be the (common) diameter of S j,n , for n ∈ N, and assume that h j → 0 as j → ∞. (For example, we might take (for n = 2) M j = {((n − 1)h j , nh j ) : n ∈ Z} × {0}.) For j ∈ N choose j > 0 with j → 0 as j → ∞, let M j denote the set of those elements of M j that have a non-empty intersection with Γ ( j ), and let Γ j denote the interior of T ∈M j T , so that M j is a convex mesh on Γ j in the sense of §A (and, in particular, is pre-convex). Then Combining Theorem 5.3 and Remark 5.4 with the density results in Lemma 3.13 we obtain the following corollary. (ii) Γ is a d-set for some n − 2 < d < n − 1 (so that V = {0}) and h j = o( µ j ) as j → ∞, for some µ > n − 1 − d.
We emphasize that in Corollary 5.6(ii) it is possible (since 0 < n − 1 − d < 1) to take µ < 1, giving convergence when h j ∼ ε j or even h j j . See the discussion around (42) and after Corollary 5.8(b) below, where this result is applied.

The BEM on fractals and prefractals arising from iterated function systems
An important IFS subclass is where each s m is a contracting similarity, i.e.
The mesh M j has N 0 elements on each component of Γ j , so N j = ν j N 0 elements in total. If h 0 is the mesh size for M 0 , then M j has mesh size h j = r j h 0 .
One can also consider meshes for which there is less than one degree of freedom (DOF) per component of Γ j . Precisely, for each j choose i = i(j) ∈ {0, . . . , j}, let τ i,1 , . . . , τ i,ν i be the components of Γ i , and consider the mesh M j on Γ j defined by ( Figure 1 shows the meshes given by (42) for j = 0, 1, 2 and 0 ≤ i ≤ j for a Cantor dust example from §6.2 below.) If i = j the mesh (42) is convex (the elements are convex sets), indeed M j coincides with the mesh given by (41) with N 0 = 1. But if i < j then the mesh M j given by (42) has only N j = ν i elements and each element is comprised of ν j−i separate components. This mesh M j is clearly not convex, if i < j, but it is pre-convex (in the sense of Appendix A) under the above assumptions on O, as captured in the following straightforward lemma. The following corollary, which follows from Theorems 4.1 and 5.3, Corollary 5.6 and Lemma 5.7, justifies convergence of the BEM when Γ j is the sequence of prefractals (39) and M j is defined by either (41) or (42), under the additional requirement that Γ ⊂ O (rather than Γ ⊂ O). We will see that this condition holds for obvious choices of O in the Cantor set and Cantor dust examples that we treat in the next sections.  (36) and that Γ ⊂ O. Then Γ is a d-set with d ∈ (0, n − 1) given by (40).
Suppose instead that M j is given by (42) with i = i(j) ∈ {0, 1, . . . , j} and let L be the diameter of Γ 0 . Then h j = r i L, and that V h j M −→ V as j → ∞ follows from Lemma 5.7 and Corollary 5.6, since, in the case n − 2 < d < n − 1, i > µj for some µ > n − 1 − d implies that h j = o( µ j ) as j → ∞.
Corollary 5.8 proves BEM convergence for the case Γ ⊂ O under rather mild mesh refinement. When d ≤ n − 2 (zero limiting solution) there is no restriction on the mesh size, in accordance with Remark 5.4. When n − 2 < d < n − 1 (nonzero limiting solution) it is possible to take µ < 1 in Corollary 5.8(ii)(b). 5 This means that BEM convergence holds with just one, or even less than one DOF per component of Γ j . In such cases, assuming that H −1/2 (Γ j ) = H −1/2 Γ j for each j, given an arbitrary , say, and then by Lemma A.1 there exists h j such that, if the mesh size on Γ j is less than h j , there exists v h j ∈ V h j such that v j − v h j H −1/2 ≤ (1 + j) −1 , and the claimed convergence follows by the triangle inequality. However, the required choice of h j depends on v j , which itself depends on v. So such an argument does not prove the existence of a single mesh refinement strategy for which V h j M −→ V . The development of a satisfactory convergence analysis for BEM on these standard prefractal sequences remains an open problem.

Examples
We now apply the theory developed above to some specific examples of fractal screens. In our first example, the Cantor set, the scattering problem is posed in R 2 (so n = 2), the screen being a subset of the one-dimensional hyperplane Γ ∞ = R × {0}. In all other examples, the scattering problem is posed in R 3 (so n = 3), the screen being a subset of the two-dimensional hyperplane Γ ∞ = R 2 × {0}. In the first three examples Γ is a compact fractal d-set (for some d < n − 1) that is the attractor (31) of some IFS of contracting similarities {s 1 , . . . , s ν }. In the remaining examples Γ is a (relatively) open subset of Γ ∞ with fractal boundary.
Interpreted in terms of DOFs, the final statement in Proposition 6.1 says that BEM convergence holds for the pre-convex mesh (42) on the thickened prefractal Γ j using (2 1− log 2 log 1/α + ) j DOFs, for arbitrary > 0. For example, for the middle third Cantor set (α = 1/3) it suffices to take 1.3 j DOFs on Γ j (note that Γ j has 2 j components).
Recalling Proposition 3.12, and comparing Propositions 6.1 and 6.2, we see that (provided the incident field doesn't vanish on the fractal screen), Cantor sets (n = 2) give rise to non-zero scattered fields for any value of α, while Cantor dusts (n = 3) give rise to non-zero scattered fields only for α > 1/4.

Sierpinski triangle
The Sierpinski triangle (or gasket) Γ is a compact subset of R 2 with empty interior, defined by (31) with ν = 3 and Since the open set condition (36) holds with O the open unit equilateral triangle with vertices (0, 0), (1, 0), (1/2, √ 3/2), Γ is a d-set with Hausdorff dimension d = log 3/ log 2. The standard prefractals Γ j are defined by (32) with Γ 0 := O (the closed unit equilaterateral triangle), so that Γ j is the (non-disjoint) union of 3 j closed equilateral triangles of side length 2 −j , and the Γ j are a decreasing nested sequence of compact sets satisfying (33). The first five prefractals are shown in Figure 3.  For the thickened prefractals, BEM convergence holds if h j = o(2 −µj ) for some µ > 2 − log 3 log 2 , in particular if h j = O(2 −j ).

Square snowflake
As our final example we consider the "square snowflake" studied in [68] (see also [40, §7.6] and the references therein). This is an open subset of R 2 with fractal boundary, constructed as the limit of a sequence of non-nested polygonal prefractals Γ j , j ∈ N 0 , the first five of which are shown in Figure 5. Each prefractal Γ j is an open polygon whose boundary is the union of N j := 4 · 8 j segments of length j := 4 −j aligned to the Cartesian axes. Let Γ 0 = (0, 1) 2 be the open unit square. For j ∈ N, ∂Γ j is constructed by replacing each horizontal edge and each vertical edge of ∂Γ j−1 respectively by the following polygonal lines composed of 8 edges each: (Note that the fourth and the fifth segments obtained are aligned; in the following however we count them as two different edges of Γ j .) Each polygonal path ∂Γ j constructed with this procedure is the boundary of a simply connected polygon Γ j of unit area, composed of 16 j squares of side length j . (See [17, §5.2] for more detail of this construction.) The resulting sequence of prefractals {Γ j } j∈N 0 is not nested: for each j ∈ N neither Γ j ⊂ Γ j−1 nor Γ j ⊃ Γ j−1 . Indeed, the two set differences Γ j \ Γ j−1 and Γ j−1 \ Γ j are composed of 4 · 8 j−1 = 2 3j−1 disjoint squares of side length j . Thus the limit set of the sequence cannot be defined simply as a union or intersection of the prefractals.
In [17] we showed how to construct inner and outer nested prefractal sequences Γ ± j such that Γ j and Γ ± j satisfy the assumptions of Proposition 4.3(iii), with the limit set Γ defined as Γ : In [17] we proved further that ∂Γ is a d-set with Hausdorff dimension d = 3/2, and that Γ is a thick domain, so that H s (Γ ) = H s Γ for all s ∈ R. Combining these facts with Theorem 3.2, Proposition 3.7, Proposition 4.3(iii) and Theorem 5.2, gives the following result. Proposition 6.5 (Square snowflake). Define the square snowflake Γ and its standard prefractals Γ j as above. Then the BVPs D(Γ ) and D(Γ j ) are well-posed, and BVP convergence holds. Furthermore, BEM convergence holds if h j → 0.

Numerical results
In this section we present numerical results validating our theory, and demonstrate the feasibility of using BEM to calculate scattering by fractal screens.
While our theoretical convergence analysis in § §5-6 is for Galerkin discretisations, the numerical results in this section were obtained using a collocation method, to make implementation as simple and flexible as possible. Our Matlab collocation code was validated against our own 2D Galerkin code for the case of the Cantor set (see Figure 8 below) and the open-source 3D Galerkin software Bempp [71] for the case of the Sierpinski triangle. In both cases, for fixed prefractal level the collocation code was found to give similar accuracy to the Galerkin codes (using the same meshes), but with a slightly lower computational cost, allowing us to reach slightly higher prefractal levels in 3D than was possible with the Galerkin code (using default Bempp settings and dense linear algebra).
All our experiments are on prefractals that are finite unions of disjoint segments (when n = 2) or finite unions of Lipschitz polygons (when n = 3). For simplicity we use uniform meshes throughout. In fact, in each experiment the elements are either congruent segments (when n = 2), or congruent squares, or congruent equilateral triangles.

Collocation method
Given a prefractal Γ j partitioned by a uniform mesh M j = {T j,1 , . . . , T j,N j } with mesh size h j (the diameter of each element of the uniform mesh), our collocation discretisation of the BIE (21) computes φ h j ∈ V h j (the N j -dimensional space of piecewise constants on M j ) by solving the equations where x l is the centre of the element T j,l . This is equivalent to approximating the testing integrals in the Galerkin equations (34) with a 1-point-per-element quadrature formula. The vector u containing the values of φ h j on each mesh element satisfies a square linear system Au = b where the matrix A and right-hand side vector b have entries A l,m = T j,m Φ(x l , x)ds(x) and b l = −g † (x l ), respectively. The integrals in the off-diagonal matrix entries are approximated with Gauss-Legendre quadrature on line segments and the tensorized version of the same rule on square elements; in both cases the number of quadrature points per element is chosen to be at least max{20h j /λ, 3} n−1 , λ = 2π/k being the wavelength. Numerical tests show that higher-order quadratures do not noticeably improve the solution accuracy for the range of parameters considered. On triangular elements we use a classical 7-point symmetric formula (as in [66, p. 415], with degree of exactness 3). The integrands of the diagonal entries A l,l have a weak singularity at the element centre x l . For line segments, A l,l is computed by dividing the segment in half and applying a high-order Gauss-Legendre quadrature on each side of the singularity. For square and equilateral triangle elements we split T j,l into 4 or 3 identical isosceles triangles respectively (with a common corner at the singularity), apply symmetry, and transform to polar coordinates, to evaluate A l,l as where L is the element side length. The integrals over θ are computed using Matlab's integral function. Since the mesh is uniform, all diagonal terms coincide and only one such computation is needed for a given value of kL. The (dense, complex, non-Hermitian) linear systems in our numerical experiments are relatively small (with fewer than 11000 DOFs) and are solved with a direct solver (Matlab's backslash).

Experiments performed
We use our BEM code to compare the numerical solutions on a sequence of prefractal screens Γ j approximating a limiting fractal screen Γ , for the examples in §6. In addition to showing domain plots of the scattered fields for different Γ j , we study the j-dependence of the norm of the numerical solution using the three norms defined in Table 1. The table also shows the marker type these norms will be represented by in all the plots. To validate our theoretical convergence results we also compute near-and far-field errors for the solution on Γ j , relative to the solution on the finest prefractal Γ j max , using these same norms. In all tests we simulate scattering problems, and the incident field is a plane wave, so that g † (x) = −e ikd·x , x ∈ Γ j , for some d ∈ S n−1 , the incident wave direction. In the convergence plots for compact screens, red continuous lines correspond to results for standard prefractals and blue dashed lines to results for thickened prefractals.

Cantor set
We first fix Γ to be the standard Cantor set, as defined in §6.1, with α = 1/3, set k = 30 (so that the acoustic wavelength λ := 2π/k ≈ 0.209 and there are roughly 5 wavelengths across Γ ), and choose the direction of the incident plane wave as d = (1/2, − √ 3/2). We make BEM computations on both the standard (open) prefractals Γ j , defined by (39) with O = (0, 1), which have 2 j components of length 3 −j (cf. §6.1), and the thickened prefractals as defined in §6.1 with δ = 1 4α − 1 2 = 1 4 , which we denote by Γ δ j , and which have 2 j components of length 3 2 3 −j . We present results for the simplest case where the BEM meshes have exactly one element per component of each prefractal, so that we are using the convex mesh (41) with N 0 = 1. Thus there are N j = 2 j elements and DOFs on the jth prefractal. For these simple meshes (Galerkin) BEM convergence is guaranteed for the thickened prefractals by Proposition 6.1. Figure 6 shows the real part and magnitude of the scattered and total fields on the box (−1, 2) × (−1.5, 1.5), computed for prefractal level j = 13, discretising Γ δ 13 with N 13 = 2 13 = 8192 elements and DOFs. The left panel in Figure 7 shows the norms, as defined in Table 1, of the (collocation) BEM solution on Γ j and Γ δ j for j = 0, . . . , 13. In all cases the norms quickly settle to an approximately constant value, suggesting that a limiting value as j → ∞ has been reached. The right panel in Figure 7 shows the near-and farfield relative errors for j = 0, . . . , 12, measured against the solutions for j = 13. Standard prefractals seem to give slightly smaller errors than thickened ones. But in both cases the relative errors decay exponentially in j, this numerical evidence of collocation BEM convergence, both for standard and thickened prefractals.  For this specific 2D problem we have also implemented a Galerkin BEM. The left panel of Figure 8 demonstrates the close agreement between our collocation solutions and the corresponding Galerkin solution (to which Proposition 6.1 applies to prove convergence in the thickened prefractal case). Taken together, since we know from Proposition 6.1 that the Galerkin solution on the thickened prefractal sequence converges to the correct limiting solution of the BIE on the Cantor set Γ , Figures 7 and 8 are persuasive numerical evidence that: i) the Galerkin solution on the standard prefractal sequence; and ii) the collocation solutions on both the standard and the thickened prefractal sequences, are all converging to the correct limiting solution for the Cantor set Γ as j → ∞. These conclusions are further supported by the right panel in Figure 8 which shows that the near and far-field relative differences between the fields computed on the standard and thickened prefractals (using collocation BEM) also decrease exponentially in j. Figure 9 shows how norms of the Cantor set solution, approximated by (collocation BEM) computations on the finest prefractal level, vary as a function of the Cantor set parameter α. Recall that α is related to the Hausdorff dimension of the Cantor set Γ by d = log 2/ log(1/α). The strength of the scattered field decreases with decreasing α (decreasing d). This is consistent with our earlier Fig. 8 Left: validation of our collocation code against a Galerkin implementation for the Cantor set problem in §7.3. Continuous red lines correspond to standard prefractals and dashed blue lines show the same quantities for thickened prefractals (see Table 1). Note that there are 6 lines in total on this graph, but the near-field and far-field relative errors are almost indistinguishable. Right: exponential decay of the relative difference between the fields scattered by the standard (Γ j ) and thickened (Γ δ j ) prefractals for the Cantor set problem in §7.3. Here φ h,δ j and u δ j,∞ denote the BEM solution on Γ δ j and the corresponding far-field pattern.

Cantor dust
We now make computations for two Cantor dusts, as defined in §6.2, with α = 1/3 and α = 1/10 respectively, setting k = 50 (so that λ ≈ 0.126 and there are roughly 11 wavelengths across the diagonal of Γ ), and choosing d = (0, 1 Similarly to §7.3 we make (collocation) BEM computations on both the standard (open) prefractals Γ j , defined by (39) with O = (0, 1) 2 (cf. §6.2), and the thickened prefractals as defined in §6.2 with δ = 1 4 , which we denote by Γ δ j . As in §7.3 we use BEM meshes with exactly one element per component of each prefractal, giving N j = 4 j DOFs in total. It follows from Proposition 6.2 that (Galerkin) BEM convergence is guaranteed for the standard prefractals for α = 1/10 (since the limiting solution is zero), and for the thickened prefractals for both α = 1/3 and α = 1/10. Figure 10 shows near-and far-field plots of the scattered field for the standard prefractal of level j = 6 with N 6 = 4 6 = 4096 DOFs.   Figure 11 shows the solution norms for prefractal levels 0 to 6 and the relative errors against the computations on the finest prefractal. From the left panels we see that for α = 1/3 the norms appear to converge to a constant positive value, while for α = 1/10 they converge (exponentially) to 0, consistent with Proposition 6.2. The superior convergence rate in the near-and far-field L 2 norms compared to the H −1/2 energy norm, visible for α = 1/10, is in line with standard superconvergence theory for functionals of a BEM solution-see e.g. [69, §4.2.5].
In the right panels of Figure 11 we observe the exponential (in j) decay of the errors of near-and far-fields against the solutions on the finest prefractal. We have also computed (but do not plot) the differences between standard and thickened prefractals in the near-and far-fields. These behave similarly to those in Figure 8 (for n = 2). These various numerical experiments, together with validations we have made (see the beginning of §7) of our collocation code against 3D Galerkin The left panel of Figure 12 shows how the magnitude of the scattered field depends on the parameter α, and thus on the Hausdorff dimension d = log 4/ log(1/α), for different prefractal levels j. Note that in this experiment we used a fixed total number N j = 4096 of DOFs on each prefractal, so that the lower order prefractal solutions are computed more accurately than they would be using our usual prescription N j = 4 j . We recall (Proposition 6.2) that in the limit j → ∞ the field should vanish for α ≤ 1/4; compare this to the right panel of Figure 9 for the Cantor set (n = 2), where the limit is non-zero for all α.
The right panel of Figure 12 shows the dependence on the wavenumber k of φ h j H −1/2 (Γ j ) , for the largest prefractal level j = 6 that approximates the fractal limit. We find that φ h 6 H −1/2 (Γ 6 ) grows with increasing k like k 0.19 for the larger values of k. In the same panel we also plot φ h Fig. 12 Left: the L 2 (S 2 ) norms of the far-field pattern for the Cantor dust, plotted against the parameter α, for standard prefractals of level j = 0, 1, . . . , 6, computed with N j = 4096 DOFs for each j. Proposition 6.2 implies that the limit for j → ∞ is 0 for α ≤ 1/4. Right: solution norms for α = 1/3 and j = 6 (in red) and j = 0 (unit square, in green) as a function of the wavenumber k.
against k (φ h 0 the numerical solution on the screen Γ 0 , i.e. the solution for a unit square screen, this computed with 10000 DOFs corresponding to more than 6 DOFs/wavelength at the highest wavenumber k = 100). Here · H −1/2 Fig. 13 Left and centre: the real part and magnitude of the field scattered by the level 8 prefractal approximation of the Sierpinski triangle, for the problem in §7.5, computed with N 8 = 3 8 = 6561 DOFs. Right: the magnitude |φ h j | of the piecewise constant BEM solution. Note that the peaks in |φ h j | (shown in yellow) are located close to the midpoints of the sides of the 9 triangular holes of size comparable to the wavelength (red segment).

Sierpinski triangle
We approximate the Sierpinski triangle with the standard prefractals Γ j described in §6.3 (to be precise we mesh Γ • j , the interior of Γ j ). We set k = 40, so that λ ≈ 0.157 and the diameter of Γ is approximately 6.4 wavelengths, and consider a downward-pointing incoming wave with d = (0, 0, −1). Figure 13 shows the near field and the magnitude of the (collocation) BEM solution φ h j for prefractal level j = 8 and N 8 = 6561 DOFs (one per component of Γ • j ). We note that |φ h j | achieves its maxima at the midpoints of the sides of the triangular holes of side length 1/8, this size comparable with the wavelength λ. Figure 14 (left panel) shows the solution norms for prefractal levels 0 to 8. In these computations we vary from our prescription in the results above of one DOF per component of the prefractal, using a larger number of elements at the lower prefractal levels to ensure that quadrature errors in the evaluation of the coefficients in the linear system are not significant with the 7-point quadrature rule that we use on triangular elements (see §7.1). Precisely, for prefractals Γ 0 to Γ 5 we use a mesh of equilateral triangles of side h j = 2 −5 . For Γ 5 to Γ 8 we use equilateral  Figure 11 for α = 1/3), we observe in the left panel rapid convergence of the norms to positive constant values. In the right panel we plot relative errors compared to the result for the finest prefractal Γ 8 . As in the right panel of Figure  11 we see convergence which is exponential in j for j ≥ 5. But the convergence is not monotonic in j for j ≤ 5 where we fix h j as we increase j.

Classical snowflakes
We now turn to examples in which the limiting screen Γ is a bounded open set with fractal boundary. In these cases, as j → ∞ the area of the prefractals must tend to the (non-zero) area of the limiting screen. Thus, in our simulations based on uniform meshes, the cost increase associated with the increasingly fine mesh required to represent the prefractal boundary exactly as j → ∞ will not be compensated by a decrease in the area to be meshed, as was the case for the Cantor sets/dusts and the Sierpinski triangle. This in turn means that our simulations based on uniform meshes and direct solvers will be more expensive than the numerical experiments reported above, limiting the prefractal level that can be attained. More efficient BEM approaches could be developed for these problems using appropriate nonuniform meshes such as those in [6,21,55], and/or fast iterative matrix-free linear solvers. But for brevity we leave such considerations to future work.
We consider first the standard Koch snowflake Γ with β = π/6 and ξ = 1/3, in the notation of §6.4. We approximate Γ with both the inner (open) and the outer (compact) prefractals Γ − j and Γ + j defined in §6.4, on which we build uniform meshes conforming to the prefractal geometries (strictly, our mesh is on (Γ + j ) • in the outer prefractal case). Figure 15 shows the scattered field for an incident plane wave with k = 61 (wavelength λ ≈ 0.103) and d = (0, 1 √ 2 , − 1 √ 2 ), at inner prefractal level 4, computed with N − 4 = 10344 DOFs. By Proposition 6.4, which relies on the fact that Γ is thick, so that H −1/2 (Γ ) = H −1/2 Γ , Galerkin BEM solutions for the inner and outer prefractals should both converge to the unique limiting solution on Γ , provided that h j → 0 as j → ∞. We Fig. 16 The relative difference, measured in the H −1/2 norm, between the BEM solutions on Γ − j in and on Γ + j out , for j in + jout = 0, 1, . . . , 7, with jout ∈ {j in − 1, j in }, for the problem described in §7.6 (green curve). We show results also for four other similar problems with different values of k and either vertically, horizontally or obliquely incident plane waves.
investigate this numerically (for the collocation BEM) in Figure 16 (though our simulations use a fixed h j on each Γ − j and a fixed mesh size also on each Γ + j , as is described in more detail below). This figure shows that the alternating inner/outer sequence of prefractal approximations has the property that the H −1/2 (R 2 ) norm of the difference between the BEM solutions on consecutive prefractals in the sequence tends to zero monotonically (and approximately exponentially) as one moves along the sequence. The figure also suggests that plane waves hitting the screen perpendicularly lead to the lowest relative difference between solutions on inner and outer prefractals, that grazing incident waves lead to the largest difference, and that the relative errors are rather insensitive to the wavenumber for the values of k investigated.
To compute an approximation to the H −1/2 (R 2 ) norm of the difference φ h j − − φ h j + , as plotted in Figure 16, we first represent both piecewise-constant fields on the same mesh. We note that the equilateral triangles of M + j are rotated by an  angle of π/2 with respect to those of M − j and are larger by a (linear) factor √ 3. We define M * j to be the smallest uniform mesh that extends M − j and covers Γ + j ; see Figure 18 for an illustration. We define two piecewise constant functions ψ ± j on this mesh. The first, ψ − j , is simply the zero-extension of φ h j − . The second, ψ + j , is defined from φ h j + as follows: noting that the centre of each element of M * j lies on an edge of M + j , we define the value of ψ + j on T * j,l ∈ M * j to be the average of the values of (the zero-extension of) φ h j + on the two triangles of (the uniform

Square snowflake
Finally, we consider the square snowflake Γ and the associated sequence of nonnested prefractals Γ j described in §6.5. We choose k = 40 and d = (0, 0, −1). Plots of the resulting near-and far-field solutions for prefractal level j = 3 are shown in Figure 19. Proposition 6.5 implies that we have (Galerkin) BEM convergence provided h j → 0 as j → ∞. But numerical validation of this convergence is hampered by the fact that the minimal number of DOFs required to discretise the jth-level prefractal with a uniform mesh of squares is 16 j . In order to simulate higher-level prefractals more sophisticated BEM implementations are needed, e.g. using fast matrix-vector multiplications, and non-uniform or non-convex meshes. This will be considered in future work.

Conclusion and open problems
In this paper we have written down in §3 well-posed BVP and BIE formulations for scattering by an acoustic screen Γ that is either fractal or has fractal boundary, refining, simplifying, and extending (in particular through application of recent results from [17]) the earlier results of [27]. Generalising the treatment in [27], we have studied in §4 (using Mosco convergence for the first time in this context) the convergence of BIE solutions on a sequence of prefractals Γ j to the solution on a limiting set Γ ; in particular we have proved in Theorem 4.1 and Proposition 4.3 new convergence results in cases where the sequence of prefractals Γ j is not monotonically convergent to Γ , and have applied these results in §4.1 to sequences of prefractals generated by general iterated function systems. But the crucial novelty of the paper has been the analysis and implementation of the BEM in § §5-7. In §4.1 we have obtained, to our knowledge, the first rigorous convergence analysis of a numerical method for scattering by a fractal object and the first proofs of convergence of BEM, applied on a sequence of prefractal sets Γ j , to the solution of the BIE on the limiting set Γ , both in cases when Γ is fractal (in particular the fixed point of an IFS of contracting similarities satisfying the standard open set condition, for example the Cantor set or dust, the Sierpinski triangle) and in cases where Γ is an open set with fractal boundary (for example the Koch or square snowflake). We have studied these specific cases as examples of applying these new convergence results in §6, and in numerical experiments in §7.
These results are, we consider, important first steps in understanding BIEs and their solutions on sets that are fractal or have fractal boundary, and the convergence to these solutions of BEM approximations on sequences of more regular prefractals. More generally, these results and methods are applicable to the convergence analysis of Galerkin methods for BIEs (and, we anticipate, other integral equations and PDEs) posed on any rough domain that cannot be discretised exactly but first needs to be approximated by a sequence of more regular sets.
For the specific screen scattering problems that have been the focus of this paper, and for large classes of related problems, there remain many intriguing and important open questions. These include: 1. What is the regularity of the BIE solution on the limiting screen Γ , and how does this depend on the fractal dimension of Γ or its boundary? In the case when Γ is fractal the density results summarised in Lemma 3.13 are one step towards addressing this question. 2. At what rate do BIE solutions on a prefractal sequence Γ j converge to the limiting BIE solution on Γ ? 3. If the BIE is additionally discretised by BEM, how does the convergence rate depend jointly on the prefractal level j and the discretisation? (The numerical simulations in §7 provide some experimental data.) 4. Extending this question, what is the optimal balance (to achieve an accurate approximation of the solution on Γ with least work) between increasing the prefractal level and mesh refinement? For example, for the family of pre-convex meshes (42), what is the optimal choice of the parameter i(j)? 5. In the case when Γ is a self-similar fractal (is the attractor of an IFS of contracting similarities), can efficient BEM solvers be built making use of the self-similarity?
We hope to address some of these questions in future work, together with exploring more efficient BEM implementations (e.g. using locally-refined meshes, fast iterative solvers for structured meshes), and the extension of the methods developed here to more general problems (e.g. curvilinear screens, different boundary conditions, elastic and electromagnetic waves, other PDEs). . We note that (ii) holds (and this will be the case in the applications that we make) if each T j is the union of a finite number of pairwise disjoint Lipschitz open sets, in particular if each T j is convex, in which case we term M a convex mesh.
In the case that Ω is a curvilinear Lipschitz polygon and each T j is a curvilinear triangle or quadrilateral, the following lemma is a standard BEM error estimate (e.g. [73,Thm. 10.4,(10.10)]), except that in the standard versions the explicit constant π s−t in (44) is replaced by an unknown constant that depends (in an unspecified way) on the domain and the shape regularity of the elements. The version we prove here, which applies to any bounded open set Ω and any pre-convex mesh M , and provides explicit constants independent of the domain and element shape, should be of some independent interest and is essential for our application to BEM on sequences of prefractals converging to a fractal limit.
Let h := max T ∈M diam(T ). Then, for −1 ≤ s ≤ 0 and 0 ≤ t ≤ 1, if u ∈ H t (Ω), then where on the left hand side we extend u − Πu from Ω to R m by zero to become an element of H 0 (Ω) ⊂ H s (Ω).