Subwavelength guided modes for acoustic waves in bubbly crystals with a line defect

The recent development of subwavelength photonic and phononic crystals shows the possibility of controlling wave propagation at deep subwavelength scales. Subwavelength bandgap phononic crystals are typically created using a periodic arrangement of subwavelength resonators, in our case small gas bubbles in a liquid. In this work, a waveguide is created by modifying the sizes of the bubbles along a line in a dilute two-dimensional bubbly crystal, thereby creating a line defect. Our aim is to prove that the line defect indeed acts as a waveguide; waves of certain frequencies will be localized to, and guided along, the line defect. The key result is an original formula for the frequencies of the defect modes. Moreover, these frequencies are numerically computed using the multipole method, which numerically illustrates our main results.


Introduction
Line defects in bandgap photonic or phononic bandgap crystals are of interest due to their possible applications in low-loss waveguides. The main mathematical problem of interest is to show that the spectrum of the defect operator has a non-zero overlap with the original bandgap. Moreover, it is also of interest to understand the nature and location of the defect spectrum. For previous works regarding line defects in bandgap crystals we refer to [11,23,24,13,14,15,16,18].
In this work, we consider a line defect in a phononic bandgap crystal comprised of gas bubbles in a liquid. The gas bubbles are known to resonate at a low frequency, called the Minnaert frequency. The corresponding wavelength is larger than the bubble by several orders of magnitude [1,30]. Based on this, it is possible to create subwavelength bandgap crystals, which operate at wavelengths much larger than the unit cell size of the microstructured material. One of the main motivations for studying subwavelength bandgap materials is to manipulate wave propagation at subwavelength scales. A second motivation is for their use in devices where conventional bandgap materials, based on Bragg scattering, would create infeasibly large devices [28,33]. Mathematical properties of bubbly phononic bandgap materials have been studied in, for example, [1,3,6,7,8,12], and subwavelength phononic bandgap materials have been experimentally realised in [25,26,27].
Wave localization due to a point defect in a bubbly bandgap material was first proven in [3]. In [2], where some additions and minor corrections to [3] were made, it is shown that the mechanism for creating localized modes using small perturbations is quite different depending on the volume fraction of the bubbles. In order to create localized modes in the dilute regime, the defect should be smaller than the surrounding bubbles, while in the non-dilute regime, the defect has to be larger. Based on this, in the case of a line defect, it is natural to expect different behaviour in these two different regimes. This suggests that different methods of analysis are needed in the two regimes. In this paper, we will mainly focus on the dilute regime, taking the radius of the bubbles sufficiently small.
If the defect size is small, i.e. if the size of the perturbed bubble is close to its original size, then the band structure of the defect problem will be a small perturbation of the band structure of the original problem [4,9]. This way, it is possible to shift the defect band upwards, and a part of the defect band will fall into the subwavelength bandgap. However, because of the curvature of the original band, it is impossible to create a defect band entirely inside the bandgap with this approach.
In order to create defect bands which are entirely located inside the subwavelength bandgap, we have to consider slightly larger perturbations. In this paper, we will show that for arbitrarily small defects, a part of the defect band will lie inside the bandgap. Moreover, we will show that for suitably large perturbation sizes, the entire defect band will fall into the bandgap, and we will explicitly quantify the size of the perturbation needed in order to achieve this. Because of this, our results are more general than previous weak localization results since we explicitly show how the defect band depends on the perturbation size.
In order to have guided waves along the line defect, the defect mode must not only be localized to the line, but also propagating along the line. In other words, we must exclude the case of standing waves in the line defect, i.e. modes which are localized in the direction of the line. As discussed in [22,23], such modes are associated with the point spectrum of the perturbed operator which appears as a flat band in the dispersion relation. Proving the absence of bound modes in phononic or photonic waveguides is a challenging problem; for example in [32] this was proven by imposing "hard-wall" Dirichlet or Neumann boundary conditions along the waveguide, while in [20] the absence of bound modes was proven in the case of a simpler Helmholtz-type operator. In this paper, we use the explicit formula for the defect band to show that it is nowhere flat, and hence does not correspond to bound modes in the direction of the line.
The paper is structured as follows. In Section 2 we discuss preliminary results on layer potentials, and outline the main results from [7]. In Section 3 we restrict to circular domains and follow the approach of [3,2] to model the line defect using the fictitious source superposition method, originally introduced in [34]. In Section 4 we prove the existence of a defect resonance frequency, and derive an asymptotic formula in terms of the density contrast in the dilute regime. Using this formula, we show that the defect modes are localized to, and guided along, the line defect. In Section 5 we compute the defect band numerically, in order to verify the formula and also illustrate the behaviour in the non-dilute regime. The paper ends with some concluding remarks in Section 6. In Appendix A, we restrict ourselves to small perturbations to derive an asymptotic formula valid in the non-dilute regime. In Appendix B we outline the fictitious source superposition method in the case of non-circular domains.

Layer potentials
Let Y 2 = [−1/2, 1/2) 2 ⊂ R 2 be the unit cell and assume that the bubble occupies a bounded and simply connected domain D ∈ Y 2 with ∂D ∈ C 1,s for some 0 < s < 1. Let Γ 0 and Γ k , k > 0 be the Green's functions of the Laplace and Helmholtz equations in dimension two, respectively, i.e., is the Hankel function of the first kind and order zero. Here, the outgoing Sommerfeld radiation condition is used for selecting the physical Helmholtz Green's function [4].
Let S k D : L 2 (∂D) → H 1 loc (R 2 ) be the single layer potential defined by Here, H 1 loc (R 2 ) denotes the space of functions that, on every compact subset of R 2 , are square integrable and have a weak first derivative that is also square integrable.
We also define the Neumann-Poincaré operator K k, * D : The following so-called jump relations of S k D on the boundary ∂D are well-known (see, for example, [4]): Here, ∂/∂ν denotes the outward normal derivative, and | ± denote the limits from outside and inside D.
In two dimensions, we have the following expansion of the Green's function for the Helmholtz equation where ln is the principal branch of the logarithm and with γ being the Euler constant. Define, for φ ∈ L 2 (∂D), Then the following expansion holds: We also introduce a quasi-periodic version of the layer potentials. For α ∈ [0, 2π) 2 , the quasi-periodic Green's function Γ α,k is defined to satisfy We define the quasi-periodic single layer potential S α,k D by It satisfies the following jump formulas: where (K −α,k D ) * is the operator given by We recall that S α,0 D :

Bubbly crystals and subwavelength bandgaps
Here we briefly review the subwavelength bandgap opening of a bubbly crystal from [7]. Assume that a single bubble occupies the region D specified in Section 2.1. We denote by ρ b and κ b the density and the bulk modulus inside the bubble, respectively. We let ρ w and κ w be the corresponding parameters outside the bubble. We introduce as the speed of sound outside and inside the bubbles, and the wavenumber outside and inside the bubbles, respectively. Here, ω corresponds to the operating frequency of the acoustic waves. Let C = ∪ n∈Z 2 (D + n) be the periodic bubbly crystal. Define, for x ∈ R 2 , where χ C is the characteristic function of C.
We assume that there is a large contrast in the density, that is, the density contrast δ satisfies Recall that under (2.3), there exists a subwavelength resonance of the bubble in free space [1].
In the following, we shall also make the assumption stated below.
Assumption 2.1. Without loss of generality, we assume that In this case we have k b = k w = ω. Assumption 2.1 only serves to simplify the expressions. The methods presented in this paper indeed apply as long as the wave speeds outside and inside the bubbles are comparable to each other.
The wave propagation problem inside the periodic crystal can be modelled as We denote by Λ 0 the set of propagating frequencies, i.e., the set of ω such that ω 2 is in the spectrum of the operator Denote by Y s = Y × R the unit strip and recall that Y 2 = [−1/2, 1/2) 2 is the unit cell of the crystal. Applying the Floquet transformation, first in x 1 -direction and then in x 2 -direction, equation (2.4) can be decomposed first as where α 1 ∈ Y * , and then as where α = (α 1 , α 2 ) ∈ Y * × Y * . We denote by Λ 0,α1 the set of ω such that ω 2 is in the spectrum of the operator implied by (2.5) and by Λ ess 0,α1 the essential part of this spectrum. It is known that (2.6) has non-trivial solutions for discrete values of ω: and we have the following band structure of propagating frequencies for the periodic bubbly crystal C: In [7], it is proved that there exists a subwavelength spectral gap opening in the band structure. Let us briefly review this result. We look for a solution v of (2.6) which has the following form: for some densities ϕ α , ψ α ∈ L 2 (∂D). Using the jump relations for the single layer potentials, one can show that (2.6) is equivalent to the boundary integral equation Since it can be shown that ω = 0 is a characteristic value for the operator-valued analytic function A(ω, 0), we can conclude the following result by the Gohberg-Sigal theory [4,19].
Theorem 2.1. [7] For α = 0 and sufficiently small δ, we have where the constant Cap D,α is given by Here, · , · stands for the standard inner product of L 2 (∂D) and χ ∂D denotes the characteristic function of ∂D.
Let ω * 1 = max α ω α 1 . The following theorem expresses the fact that a subwavelength bandgap opens in the band structure of the bubbly crystal.
3 Integral representation for bubbly crystals with a defect

Formulation of the line defect problem
In the following, we will consider the case when all the bubbles are circular disks. This gives a convenient presentation, and makes the problem similar to the point defect problem studied in [3,2]. In Appendix B, we will outline the analysis in the case of non-circular bubbles. Consider a perturbed crystal, where all the disks along the x 1 -axis are replaced by defect disks of radius R d with 0 < R d < R. Denote the centre defect disk by D d and let The wave propagation problem inside the periodic crystal can be modelled as We denote by Λ d the set of propagating frequencies in the line defect crystal, i.e. the set of ω such that ω 2 is in the spectrum of the operator Since the defect crystal is periodic in the x 1 -direction, we can use the Floquet transformation to decompose (3.1) as where α 1 ∈ Y * and Y s again denotes the strip Y s = [−1/2, 1/2) × R. We will denote by Λ d,α1 the set of ω such that ω 2 is in the spectrum of the operator implied by (3.2) and by Λ ess d,α1 the corresponding essential part of the spectrum.
In the strip Y s , the perturbations ρ d −ρ and κ d −κ have compact support. Since the essential spectrum is stable under compact perturbations [17,31], it can be shown that the essential spectra Λ ess 0,α1 and Λ ess d,α1 coincide.
In this paper, we want to show that introducing the line defect creates a defect band ω ε (α 1 ) / ∈ Λ 0,α1 . Moreover, we want to show that ε can be chosen such that ω ε (α 1 ) / ∈ Λ 0 for all α 1 ∈ Y * , which means that any Bloch mode is localized to the line defect. We also want to show that ω ε (α 1 ) is not contained in the pure point part of Λ d,α1 , which means that there are no bound modes in the defect direction.

Effective sources for the defect
Here we describe an effective sources approach to the solution of (3.2) in the strip. The idea is to model the defect bubble D d as an unperturbed bubble D with additional fictitious monopole and dipole sources f and g. This method was originally introduced in [34] and then it was applied in [3,2] for a point defect in a bubbly crystal.
Let us consider the following problem: where f and g are the source terms and δ m,n is the Kronecker delta function. Note that the sources are present only on the boundary of the central bubble D.
We denote the solution to the original problem (3.2) by u and the effective source solution (3.3) by u. We want to find appropriate conditions on f and g in order to achieve Then u can be recovered by extending u to the whole region including D \ D d with boundary conditions on ∂D and ∂D d . The conditions for the effective sources f and g, which are necessary in order to correctly model the defect, will be characterized in the next subsection.

Characterization of the effective sources
Here we clarify the relation between the effective source pair (f, g) and the layer density pair (ϕ, ψ) defined in equation (3.5) below. First, we observe that away from the central unit cell Y 2 , the equations (3.2) and (3.3) satisfy the same geometric and quasi-periodic conditions. Thus, in order for (3.4) to hold, it is sufficient for u and u to coincide inside the central unit cell Y 2 .
Inside Y 2 , the solution u can be represented as in D, (3.5) for some pair (ϕ, ψ) ∈ L 2 (∂D) 2 , with H satisfying the homogeneous equation (∆ + k 2 w )H = 0 in Y 2 . In (3.5), the local properties of u around ∂D are given by the single-layer potentials, while H can be chosen to make u satisfy the quasi-periodic condition. From the jump conditions given in Section 2.1, the pair (ϕ, ψ) satisfies Similarly, inside Y 2 , the solution u can be represented as Now, having the two solutions coincide inside (Y 2 \ D) ∪ D d is equivalent to the conditions and Assuming D is a disk, the above equations were solved in [3,2], and we state the results in Proposition 3.1 below. First, we introduce some notation. Since D and D d are circular disks, we can use a Fourier basis for functions in L 2 (∂D) or L 2 (∂D d ). For n ∈ Z, define the subspace V n of L 2 (∂D) as V n := span{e imθ }. Then define the subspace V mn of L 2 (∂D) 2 as Similarly, let V d mn be the subspace of L 2 (∂D d ) 2 with the same Fourier basis. Then it can be shown that the operator A D in (3.6) has the following matrix representation as an operator from V mn to V m ′ n ′ : Similarly, the operator A D d in (3.7) is represented as follows: In [3,2], the following proposition was shown.
Proposition 3.1. The density pair (ϕ, ψ) and the effective sources (f, g) satisfy the following relation where the operators P 1 :

Floquet transform of the solution
In view of Proposition 3.1, we can identify the solutions u and u. In this section, we derive an integral equation for the effective source problem (3.3). This problem is already quasi-periodically reduced in the x 1 -direction, with quasi-periodicity α 1 . For some quasi-periodicity α 2 ∈ Y * , we set α = (α 1 , α 2 ) and apply the Floquet transform to the solution u in the x 2 -direction as follows: The transformed solution u α satisfies The solution u α is α-quasi-periodic in the two-dimensional cell Y 2 , and can be represented using quasiperiodic layer potentials as where, similarly as in equation (2.7), the pair (ϕ α , ψ α ) ∈ L 2 (∂D) 2 is the solution to Since the operator A α is invertible for small enough δ and for ω inside the bandgap [7], we have The solution u to problem (3.3) can be recovered by the inversion formula as Now, by the same arguments as those in [3,2], we obtain the following proposition.

The integral equation for the layer densities
Here we state the integral equation for the layer density pair (φ, ψ). The following result is an immediate consequence of Propositions 3.1 and 3.2.
The expression of this integral equation resembles the one for a point defect found in [3,2]. However, this similarity is not obvious, and can be seen as a consequence of the cancellation of H in Proposition 3.1.
The significance of Proposition 3.3 is as follows. If we can show that there is a characteristic value ω = ω ε of M ε,δ,α1 inside the bandgap, i.e. if there is a non-trivial pair (φ, ψ) such that M ε,δ,α1 (ω ε ) φ ψ = 0, then ω ε is a resonance frequency for the defect mode.

Subwavelength guided modes in the defect
Here, we will prove the existence of a resonance frequency ω = ω ε (α 1 ) inside the bandgap of the unperturbed crystal at α 1 . We will give an asymptotic formula for M ε,δ,α1 in terms of δ in the dilute regime. Moreover, we will show that the defect band is not contained in the pure point spectrum of the defect operator, and for perturbation sizes ε larger than some critical ε 0 , the entire defect band is located in the bandgap region of the original operator.

Define the operators
The motivation for defining these operators is given in Lemmas 4.1 and 4.2. Introducing these operators enables the explicit computation of (A α ) −1 . We will compute the asymptotic expansion of these operators for small ω and δ.
Proof of (ii). Using the definition of A 1 , and the expression for A 0 , we can compute Recall the asymptotic expression of ω α given in Theorem 2.1: We then find that .

Lemma 4.2.
For ω = ω α , and as ω → 0 and δ → 0, we have Proof. We have already established the invertibility of A 0 and A 1 . Using this fact, we have Because Q * α χ ∂D = 0, we have that We then have that where the last step follows using the Neumann series.
Next, we compute the operator (A ε D − A D ). Using Proposition 3.1 and equation (3.10), we have .
Consequently, the operator (A ε D − A D ) is given by

Introduce the notation
(4.6) Using asymptotic expansions of the Bessel function J n (z) and the Hankel function H n (z), for small z, straightforward computations show that Moreover, we have We are now ready to compute the full operator M ε,δ,α1 .

Defect resonance frequency in the dilute regime
The following theorem is the main result of this paper. Again, we say a frequency ω is in the subwavelength regime if ω = O( √ δ).
Theorem 4.1. For δ and R small enough, there is a unique characteristic value ω ε (α 1 ) of M ε,δ,α1 (ω) such that ω ε (α 1 ) = Λ 0,α1 and ω ε (α 1 ) is in the subwavelength regime. Moreover, whereω is the root of the following equation: Proof. We seek the characteristic values of the operator I + M 0 . We consider the dilute regime, i.e. where R is small. As shown in [7], in this case we have where R α (x) = Γ α,0 (x) − Γ 0 (x). In particular, We will compute M 0 in the Fourier basis. It is known that [4]  In total, we have on the subspace V 0 , Moreover, if n = 0, then In summary, the operator I + M 0 can be written as where the limiting operatorM (ω) is a diagonal operator in the Fourier basis, with non-zero diagonal entries for n = 0. We conclude that ω =ω is a characteristic value forM (ω) if and only if one of the diagonal entries vanishes at ω =ω, i.e. if (4.11) Next, we show that equation (4.11) has a zeroω / ∈ Λ 0,α1 satisfyingω = O( √ δ). Introduce the notation For a fixed α 1 ∈ Y * , define ω * = ω (α1,π) , which is the edge of the first band in Λ 0,α1 . Observe that 1/I(ω, α 1 ) is monotonically increasing in ω, and where ω 2 0 is the average In the dilute regime, we can compute so as ω → ω * , the right-hand side of equation (4.12) tends to Since R d < R, the leading-order term is negative. On the other hand, as ω → ∞, the right-hand side of (4.12) tends to ∞. Since the right-hand side of equation (4.12) is monotonically increasing, this equation has a unique zero ω =ω. It can be verified that this zero has multiplicity one. Moreover,ω satisfieŝ ω = O( √ δ). Now, we turn to the full operator I + M 0 (ω). Since I + M 0 (ω) =M (ω)+ O R 2 + ω , by the Gohberg-Sigal theory [4,9,19], close toω there is a unique characteristic value ω ε of the operator I + M 0 (ω), This concludes the proof.

Remark 4.2.
In the case of R d > R, i.e. larger defect bubbles, similar arguments show that any subwavelength frequency ω ε (α 1 ) ∈ Λ d,α1 \ Λ 0,α1 satisfies equation (4.9) in the dilute regime. However, it is easily verified that this equation has no solutionsω > ω * in the case R d > R. The conclusion is that we must reduce the size of the defect bubbles in order to create subwavelength guided modes in the dilute regime.

Absence of bound modes in the line defect direction
In this section, we will show that the defect band is not contained in the pure point spectrum of the defect crystal.
Proof. From the spectral form of the Green's function [4]: it can be easily shown that By symmetry of the summation, we find that  For δ and R small enough, and for α 1 = 0, π, the characteristic value ω ε = ω ε (α 1 ) satisfies ∂ω ε ∂α 1 = 0.
Proof. To simplify the computations, we introduce the following notation: Then equation (4.9) reads Denote by x ′ = ∂x ∂α1 and y ′ = ∂y ∂α1 , then we have First, we show that A = 0 which implies that the zeros of x ′ coincides with the zeros of B. We have Next, we show that the leading order of B vanishes exactly at the points α 1 = 0 and α = π. Using equations (4.5) and (4.10), we have Since R α = Γ α,0 − Γ 0 , using Lemma 4.3 we conclude that for δ and R small enough and α 1 = 0, π, y ′ is non-zero for any α 2 . Hence B is non-zero, which concludes the proof.
Proposition 4.2 shows that the defect dispersion relation is not flat, except for local extrema at α 1 = 0 and α 1 = π. Thus, the defect band is not in the pure point spectrum of the defect operator, and corresponding Bloch modes are not bounded in the line defect direction.

Bandgap located defect bands
In this section, we will demonstrate that it is possible to position the entire defect band in the bandgap region with a suitable choice of ε. Recall that ε = R d − R. As before, let is attained at α 1 = α 0 with α 0 → 0 as δ → 0.

Discretization of the operator
The operator M ε,δ,α1 (ω) was approximated as a matrix M (ω) using the truncated Fourier basis e −iN θ , e −i(N −1)θ , . . . , e iN θ . We refer to [3,2] for the details of the discretization. The integral over Y * in (3.12) was approximated using the trapezoidal rule with 100 discretization points. The characteristic value problem for M ε,δ,α1 (ω) was formulated as the root-finding problem det M (ω) = 0 and solved using Muller's method [4].

Evaluation of the asymptotic formula
The integral over Y * in equation (4.9) was approximated using the trapezoidal rule with 100 discretization points. Again, the equation was numerically solved using Muller's method. Figure 2 shows the unperturbed band structure and the defect band for α 1 over the Brillouin zone [0, 2π].

Dilute regime
The material parameters were chosen as κ b = ρ b = 1, κ w = ρ w = 5000, R = 0.05 and ε = −0.2R. It can be seen that the entire defect band is located inside the deep subwavelength regime of the bandgap. Moreover, the defect frequencies computed using the asymptotic formula agree well with the values computed by discretizing the operator M ε,δ,α1 . Also, we see that the defect band is not flat. In summary, these results show that the defect crystal supports guided modes in the subwavelength regime, localized to the line defect.

Computation of ε 0
In this section, we numerically compute the critical perturbation size, where the entire defect band is located in the bandgap. The critical perturbation size was computed in two ways: by solving equation (4.13) for the leading order term, and by solving the root-finding problem ω ε0 (0) = ω * where ω ε was computed by discretizing the operator M ε,δ,α1 . Figure 3 shows ε 0 for different R in the dilute regime. The material parameters were chosen as κ b = ρ b = 1 and κ w = ρ w = 10000. The values obtained from the asymptotic formula and by discretizing the operator agree, with a smaller radius R giving a smaller error. Quantitatively, for R in this regime, we require a decrease of the bubble size by around 14% to 26% in order that the defect band be located inside the bandgap.

Non-dilute regime
Here we compute the defect band in the non-dilute regime, in both cases ε < 0 and ε > 0, corresponding to smaller and larger defect bubbles, respectively. Theorem A.1 in Appendix A shows that there is a defect frequency ω ε in the bandgap for small ε > 0 but not for small ε < 0. Figure 4 shows the band structure in the non-dilute case with ε > 0. The material parameters were chosen as κ b = ρ b = 1, κ w = ρ w = 5000, R = 0.4 and ε = 0.45R. As expected from Theorem A.1, there is a defect band above the first band of the unperturbed crystal. Moreover, it is possible to position the entire band inside the bandgap.   Figure 5 shows the band structure in the non-dilute case with ε = −0.6R. The material parameters were chosen as κ b = ρ b = 1, κ w = ρ w = 5000, R = 0.4 and ε = 0.45R. In this case a defect band is present inside the bandgap. Here ε is quite large, in contrast to Theorem A.1 which is only valid for small ε.

Concluding remarks
In this paper, we have for the first time proved the possibility of creating subwavelength guided waves localized to a line defect in a bubbly phononic crystal. We have shown that introducing a defect line, by shrinking the bubbles along the line, creates a defect frequency band inside the bandgap of the original crystal. An arbitrarily small perturbation will create a non-zero overlap between the defect band and the bandgap, and we have explicitly quantified the required defect size in order to position the entire defect band inside the bandgap. Moreover, we have shown for the first time that the defect band is not contained in the pure point spectrum of the perturbed operator. This shows that we can create truly guided modes, which are not localized in the direction of the defect. In the future, we plan to study more sophisticated waveguides, with bends and junctions. Moreover, we also plan to study waveguides in phononic subwavelength bandgap crystals with non-trivial topology, rigorously proving the existence of topologically protected subwavelength states in bubbly crystals.
A The resonance frequency of the defect mode for small perturbations Here we derive a formula for the resonance frequency of the defect mode in the case of small ε, following the approach of [3,2]. The strength of this approach is that it is valid in both the dilute and non-dilute regimes. We begin by reformulating the integral equation (3.12) in terms of the effective sources (f, g) instead of the layer densities (φ, ψ). The following proposition is a restatement of Proposition 3.3.
In this section, we derive an expression for the characteristic value ω ε of M ε,δ,α1 (ω) located slightly above ω α for both the dilute and non-dilute regimes.
Let us first analyse the operator Y * (A α ) −1 dα. Since ω α is a simple pole of the mapping ω → A α (ω, ) −1 in a neighbourhood of ω α , according to [4], we can write where the operator-valued function R α (ω) is holomorphic in a neighbourhood of ω α , and the operator L α maps L 2 (∂D) 2 onto ker A α (ω α , δ). Let us write where * denotes the adjoint operator. Then, as in [4,5], it can be shown that where again · , · stands for the standard inner product of L 2 (∂D) 2 . Hence the operator M ε,δ,α1 can be decomposed as Note that the third term in the right-hand side is holomorphic with respect to ω. Denote by α * = (α 1 , π) and ω * (α 1 ) = ω (α1,π) . Using similar arguments as in [10] and the fact that each bubble is a circular disk, we can prove the following result on the shape of the dispersion relation close to α * .
Lemma A.2. The following results hold: which is positive for δ small enough.

Combining equation (A.3) and Lemma A.2, we obtain the following result
Theorem A.1. Assume that δ is small enough and the pair (R, ε) satisfies one of the two assumptions (i) R small enough and ε < 0 small enough in magnitude (Dilute regime), (ii) R close enough to 1/2 and ε > 0 small enough (Non-dilute regime).
Remark A.1. It is easily verified that equation (4.9) evaluated for small ε coincides with equation (A.4) evaluated for small R.

B Characterization of the effective sources for non-circular bubbles
Let now D be a general simply connected domain with ∂D ∈ C 1 . In this section, we will restrict to the case of small size perturbations ε < 0. Define the defect bubble D d ∈ D as the domain with boundary where ν x is the outward unit normal of ∂D at x ∈ ∂D. We will need some results given in [5]. First, we introduce some notation. Define the mapping p : ∂D → ∂D d , p(x) = x + εν x . Let x, y ∈ ∂D and let x = p(x) ∈ ∂D d and y = p(y) ∈ ∂D d . Define q : L 2 (∂D) → L 2 (∂D d ), q(φ)( x) = φ(p −1 ( x)), and for a surface density φ on ∂D, define φ = q(φ) on ∂D d . We also define the signed curvature τ = τ (x), x ∈ ∂D in the following way. Let x = x(t) be a parametrization of ∂D by arc length. Then define τ by Observe that τ is independent of the orientation of ∂D. The following results are given in [5], but adjusted to the case where ε < 0.
Proposition B.1. Let k > 0. Let φ ∈ L 2 (∂D) and let x, y, x, y, φ be as above. Then Proposition B.2. Let φ ∈ L 2 (∂D) and let x, y, x, y, φ be as above. Then where K k 1 is given by Here ∂ 2 ∂T 2 denotes the second tangential derivative, which is independent of the orientation of ∂D. We also state the following result which is given, for example, in [29].
As in Section 4, we consider the defect problem (3.1), modelled by the fictitious sources as in equation (3.3). Observe that Proposition 3.2 is valid even for the case of non-circular bubbles. To derive the analogue of Proposition 3.1, we again study equations (3.8) and (3.9), i.e., Since ω is in the subwavelength regime, k b is not a Dirichlet eigenvalue. Together with the uniqueness of the exterior Dirichlet problem, we conclude that it is sufficient to consider these equations on the boundaries. Using the notation from above, this means Using the expansions (B.1), (B.2) and (B.3), we find with q defined as above. From this we find that where Q is the bijection Q : L 2 (∂D) 2 → L 2 (∂D d ) 2 , Q = (q, q) and P 1 : L 2 (∂D) 2 → L 2 (∂D d ) 2 .
Using Using Taylor expansion, we have that We use the Laplacian in the curvilinear coordinates defined by T x , ν x for x ∈ ∂D d , It is easily verified that the curvatures on the two boundaries satisfy τ ( x) = τ (x) + O(ε).
Proposition B.4. The density pair (ϕ, ψ) and the effective sources (f, g) satisfy the following relation where the operator A ε D satisfies