---
abstract: 'The ability to control wave propagation is of fundamental interest in many areas of physics. Photonic and phononic crystals have proved very useful for this purpose but, because they are based on Bragg interference, these artificial media require structures with large dimensions. In \[Ammari et al., Subwavelength phononic bandgap opening in bubbly media, J. Diff. Eq., 263 (2017), 5610–5629\], it has been proved that a subwavelength bandgap opening occurs in bubbly phononic crystals. To demonstrate the opening of a subwavelength phononic bandgap, a periodic arrangement of bubbles is considered and their subwavelength Minnaert resonance is exploited. In this paper, this subwavelength bandgap is used to demonstrate cavities, very similar to those obtained in photonic and phononic crystals, albeit of deeply subwavelength dimensions. The key idea is to perturb the size of a single bubble inside the crystal, thus creating a defect. The goal is then to analytically and numerically show that this crystal has a localized eigenmode close to the defect bubble.'
author:
- 'Habib Ammari[^1]'
- Brian Fitzpatrick
- 'Erik Orvehed Hiltunen[^2]'
- Sanghyeon Yu
bibliography:
- 'defect\_final.bib'
title: Subwavelength localized modes for acoustic waves in bubbly crystals with a defect
---

2 2

2[35R30, 35C20.]{}

Introduction
============

It is well-known in solid state physics that the periodicity of atoms composing crystals is responsible for the existence of both conducting bands and bandgaps for electrons. This property is a consequence of the Floquet-Bloch theory applied to the wave-function of electrons. Similarly, electromagnetic and elastic waves propagating in periodic media are subject to the same formalism giving rise to the existence of ranges of frequencies for which no propagation is allowed, so called bandgaps. Such materials are known as photonic and phononic crystals. At an interface with free space, these types of crystals act as mirrors for incoming waves, termed Bragg mirrors. Similarly, a point defect can be created by locally modifying the properties of a crystal. This results in a cavity if the mode supported by the defect is resonant within the bandgap. Physically, the underlying mechanism is Bragg interference. In these periodic media, the Bragg condition typically occurs when the period of the medium scales with the wavelength [@arma; @lecturenotes; @soussi; @figotinkuchment; @kuchment2; @hempel; @Lipton; @Lipton2]. As a consequence, photonic and phononic crystals are typically structured with a period corresponding to half the operating wavelength. This constrains the range of applications, specifically in low-frequency regimes where the wavelength is large [@pnas; @nature].

Based on the phenomenon of subwavelength resonance, a class of phononic crystals that exhibit bandgaps with deep subwavelength spatial scales have been fabricated [@phononic1]. In [@bandgap], the opening of a subwavelength phononic bandgap in bubbly crystals has been proved. This subwavelength bandgap is mainly due to the cell resonance of the bubbles in the quasi-static regime and is quite different from the usual bandgaps in photonic/phononic crystals where the gap opens at wavelengths which, as mentioned previously, are comparable to the period of the structure. In [@highfrequency], it has been further proved that the first Bloch eigenvalue achieves its maximum at the corner of the Brillouin zone. Moreover, by computing the asymptotic behaviour of the Bloch eigenfunctions in the periodic structure near that critical frequency, it has been demonstrated that these eigenfunctions can be decomposed into two parts: one which is slowly varying and satisfies a homogenized equation, while the other is periodic across each elementary crystal cell and is varying more rapidly. The asymptotic analysis of wave fields near the critical frequency where a subwavelength bandgap opens rather than the zero frequency has been performed. This rigorously justifies the observed superfocusing of acoustic waves in bubbly crystals near and below the maximum of the first Bloch eigenvalue and confirms the bandgap opening near and above this critical frequency. We refer the reader to [@superfocusing; @Ammari2015; @Ammari2015_a] for the mathematical analysis of the superfocusing phenomenon in resonant media.

Bubbly media is a natural model for the control of wave propagation at the deep subwavelength scale because of the simplicity of the constituent resonant structure, the air bubbles. It is well-known that a single bubble in water possesses a subwavelength resonance which is called the Minnaert resonance [@first; @Minnaert]. By the subwavelength resonance, we mean that the resonator (in our case, the bubble) is of size smaller than the the wavelength. This resonance is due to the high contrast in density between the bubble and the surrounding medium and it makes the air bubble an ideal subwavelength resonator (the bubble can be two order of magnitude smaller than the wavelength at the resonant frequency).

In phononic crystals, a point defect can be created by locally removing a scatterer. This results in a small cavity because a resonant mode is created by the defect within the bandgap [@weinstein2; @weinstein3; @figotin1; @figotin2; @weinstein1; @klein; @lin1; @lin2]. Following this concept, many components have been demonstrated based on periodic media such as waveguides using line defects [@santosa; @santosa2]. However, because of their wavelength scale period, phononic crystals result in relatively large devices. This seriously constrains the range of applications, especially in the low frequency regime where the wavelength is large.

If we remove one bubble inside the bubbly crystal, we cannot create a defect mode. The defect created in this fashion is actually too small to support a resonant mode, while in a phononic crystal removing one scatterer allows for the existence of a stationary defect mode since the typical scale of such a defect is the wavelength. This illustrates a strong difference between Bragg bandgaps and subwavelength bandgaps in bubbly crystals.

In order to tackle this issue, we have to physically introduce a resonant defect inside the crystal of subwavelength resonators, and this is achieved by simply detuning one resonator with respect to the rest of the medium. In the case of the bubbly medium, we prove in this paper that by perturbing the radius of one bubble we create a detuned resonator with a resonance frequency that is upward shifted, thus falling within the subwavelength bandgap. Moreover, we will show that the way to shift the frequency upwards depends on the crystal: in the dilute regime we have to decrease the defect bubble size while in the non-dilute regime we have to increase the size.

The aim of this paper is to prove the existence of this defect mode. Through the application of layer potential techniques, Floquet-Bloch theory, and Gohberg-Sigal theory we derive an original formula for the defect mode frequency, along with proving the existence of a subwavelength localized mode. Our results are complemented by several numerical examples which serve to validate them in two dimensions. Our results formally explain the experimental observations reported in [@experiment2013; @experiment] in the case of Helmholtz resonators. They lay the mathematical foundation for the analysis of wave propagation control the deeply subwavelength scale. Subwavelength cavities have a high quality factor and a low mode volume. These two effects are typically associated with the enhancement of the emission rate of an emitter or the so-called Purcell factor [@purcell].

The paper is organized as follows. In Section \[sec-1\] we formulate the spectral problem for a bubble phononic crystal and introduce some basic results regarding the quasi-periodic Green’s function, stability of the essential spectrum, and Floquet-Bloch theory. In Section \[sec-2\] we use the fictitious source superposition method introduced in [@Wilcox] for modelling the defect and characterize the fictitious sources as the solution of some system of integral equations. In Section \[sec-3\], we prove existence of a localized defect mode and derive an asymptotic formula for the resonant frequency created inside the subwavelength bandgap by perturbing the size of a single bubble in terms of the difference between its radius and the radius of the original bubbles. In Section \[sec-4\] we perform numerical simulations to illustrate the main findings of this paper. We make use of the multipole expansion method to compute the defect mode inside the subwavelength bandgap. The paper ends with some concluding remarks.

Preliminaries {#sec-1}
=============

Layer potentials
----------------

Let $\Gamma^0$ and $\Gamma^k,k>0$ be the fundamental solution of the Laplace and Helmholtz equations in dimension two, respectively, [*i.e.*]{}, $$\begin{cases}
{\displaystyle}\Gamma^k(x,y) = -\frac{i}{4}H_0^{(1)}(k|x-y|), \ & k>0, \\
{\noalign{\smallskip}}{\displaystyle}\Gamma^0(x,y) = \frac{1}{2\pi}\ln|x-y|,
\end{cases}$$ where $H_0^{(1)}$ is the Hankel function of the first kind and order zero.

Let ${\mathcal{S}}_{D}^k: L^2(\partial D) \rightarrow H_{\textrm{loc}}^1({\mathbb{R}}^2)$ be the single layer potential defined by $${\mathcal{S}}_D^k[{\varphi}](x) = \int_{\partial D} \Gamma^k(x,y){\varphi}(y) {\: \mathrm{d}}\sigma(y), \quad x \in {\mathbb{R}}^2.$$ We also define the Neumann-Poincaré operator $\left({\mathcal{K}}_D^k\right)^*: L^2(\partial D) \rightarrow L^2(\partial D)$ by $$\left({\mathcal{K}}_D^k\right)^* = \int_{\partial D} \frac{\partial }{\partial \nu_x}\Gamma^k(x,y) {\varphi}(y) {\: \mathrm{d}}\sigma(y), \quad x \in \partial D.$$ In the case when $k=0$, we will omit the superscripts and write ${\mathcal{S}}_D$ and ${\mathcal{K}}_D^*$, respectively. The following so-called jump relations of ${\mathcal{S}}_D^k$ on the boundary $\partial D$ are well-known (see, for instance, [@lecturenotes]): $${\mathcal{S}}_D^k[{\varphi}]\big|_+ = {\mathcal{S}}_D^k[{\varphi}]\big|_-,$$ and $$\frac{\partial }{\partial \nu}{\mathcal{S}}_D^k[{\varphi}] \bigg|_{\pm} =  \left(\pm\frac{1}{2} I + \left({\mathcal{K}}_D^k\right)^*\right) [{\varphi}].$$ Here, $\partial/\partial \nu$ denotes the outward normal derivative, and $|_\pm$ denotes the limits from outside and inside $D$. In two dimensions, the fundamental solution of the free-space Helmholtz equation has a logarithmic singularity given by [@lecturenotes] $$\label{eq:hankel}
-\frac{i}{4}H_0(k|x-y|) = \frac{1}{2\pi} \ln |x-y| + \eta_k + \sum_{j=1}^\infty\left( b_j \ln(k|x-y|) + c_j \right) (k|x-y| )^{2j},$$ where $\ln$ is the principal branch of the logarithm and $$\eta_k = \frac{1}{2\pi}(\ln k+\gamma-\ln 2)-\frac{i}{4}, \quad b_j=\frac{(-1)^j}{2\pi}\frac{1}{2^{2j}(j!)^2}, \quad c_j=b_j\left( \gamma - \ln 2 - \frac{i\pi}{2} - \sum_{n=1}^j \frac{1}{n} \right),$$ and $\gamma$ is the Euler constant. Define $$\hat{{\mathcal{S}}}_D^k[{\varphi}](x) = {\mathcal{S}}_D[{\varphi}](x) + \eta_k\int_{\partial D} {\varphi}{\: \mathrm{d}}\sigma.$$ Then the following expansion holds: $$\label{eq:Sexpansion}
{\mathcal{S}}_D^k =  \hat{{\mathcal{S}}}_{D}^k + {\mathcal{O}}(k^2 \ln k).$$

We also introduce a quasi-periodic version of the layer potentials. Let $Y$ be the unit cell $[-1/2,1/2)^2$ in $\mathbb{R}^2$. For $\alpha\in [0,2\pi)^2$, the function $\Gamma^{\alpha, k}$ is defined to satisfy $$(\Delta_x + k^2) \Gamma^{\alpha, k} (x,y) = \sum_{n\in \mathbb{R}^2} \delta(x-y-n) e^{i n\cdot \alpha},$$ where $\delta$ is the Dirac delta function and $\Gamma^{\alpha, k} $ is $\alpha$-quasi-periodic, [i.e.]{}, $e^{- i \alpha\cdot x} \Gamma^{\alpha, k}(x,y)$ is periodic in $x$ with respect to $Y$.

We define the quasi-periodic single layer potential $\mathcal{S}_D^{\alpha,k}$ by $$\mathcal{S}_D^{\alpha,k}[{\varphi}](x) = \int_{\partial D} \Gamma^{\alpha,k} (x,y) {\varphi}(y) d\sigma(y),\quad x\in \mathbb{R}^2.$$ It satisfies the following jump formulas: $${\mathcal{S}}_D^{\alpha,k}[{\varphi}]\big|_+ = {\mathcal{S}}_D^{\alpha,k}[{\varphi}]\big|_-,$$ and $$\frac{{\partial}}{{\partial}\nu} \Big|_{\pm} \mathcal{S}_D^{\alpha,k}[{\varphi}] = \left( \pm \frac{1}{2} I +( \mathcal{K}_D^{-\alpha,k} )^*\right)[{\varphi}]\quad \mbox{on}~ {\partial}D,$$ where $(\mathcal{K}_D^{-\alpha,k})^*$ is the operator given by $$(\mathcal{K}_D^{-\alpha, k} )^*[{\varphi}](x)= \mbox{p.v.} \int_{{\partial}D} \frac{{\partial}}{{\partial}\nu_x} \Gamma^{\alpha,k}(y,y) {\varphi}(y) d\sigma(y).$$ We remark that it is known that $\mathcal{S}_D^{\alpha,0} : L^2({\partial}D) \rightarrow H^1({\partial}D)$ is invertible for $\alpha \ne 0$ [@lecturenotes].

Floquet transform
-----------------

A function $f(x)$ is said to be $\alpha$-quasi-periodic in the variable $x\in {\mathbb{R}}^2$ if $e^{-i\alpha\cdot x}f(x)$ is periodic. Given a function $f\in L^2({\mathbb{R}}^2)$, the Floquet transform is defined as $$\label{eq:floquet}
{\mathcal{F}}[f](x,\alpha) = \sum_{m\in {\mathbb{Z}}^2} f(x-m) e^{i\alpha \cdot m},$$ which is $\alpha$-quasi-periodic in $x$ and periodic in $\alpha$. Let $Y = [-1/2,1/2)^2$ be the unit cell and $BZ := {\mathbb{R}}^2 / 2\pi {\mathbb{Z}}^2 \simeq [0,2\pi)^2$ be the so-called first Brillouin zone. The Floquet transform is an invertible map ${\mathcal{F}}:L^2({\mathbb{R}}^2) \rightarrow L^2(Y\times BZ)$, with inverse (see, for instance, [@kuchment; @lecturenotes]) $${\mathcal{F}}^{-1}[g](x) = \frac{1}{(2\pi)^2}\int_{BZ} g(x,\alpha) {\: \mathrm{d}}\alpha.$$

Bubbly crystals and subwavelength bandgaps {#subsec:bandgap}
------------------------------------------

Here we briefly review the subwavelength bandgap opening of a bubbly crystal from [@bandgap].

Assume that a single bubble occupies $D$, which is a disk of radius $R<1/2$ centred at the origin. We denote by $\rho_b$ and $\kappa_b$ the density and the bulk modulus of the air inside the bubble, respectively. We let $\rho_w$ and $\kappa_w$ be the corresponding parameters for the water. We introduce $$v_w = \sqrt{\frac{\kappa_w}{\rho_w}}, \quad v_b = \sqrt{\frac{\kappa_b}{\rho_b}}, \quad k_w= \frac{\omega}{v_w} \quad \text{and} \quad k_b= \frac{\omega}{v_b}$$ to be the speed of sound outside and inside the bubbles, and the wavenumber outside and inside the bubbles, respectively. $\omega$ corresponds to the operating frequency of acoustic waves. We also introduce two dimensionless contrast parameters $$\delta = \frac{\rho_b}{\rho_w} \quad \text{and} \quad \tau= \frac{k_b}{k_w}= \frac{v_w}{v_b} =\sqrt{\frac{\rho_b \kappa_w}{\rho_w \kappa_b}}.$$ We assume that the wave speeds outside and inside the bubbles are comparable to each other and that there is a large contrast in the bulk moduli, that is, $$\delta \ll 1,\quad \tau= O(1).$$ In this paper, for the sake of simplicity of presentation, we shall assume $v_w=v_b=1$ .

Let ${\mathcal{C}}= \cup_{n\in{\mathbb{Z}}^2}(D+n)$ be the periodic bubbly crystal. Consider the following quasi-periodic scattering problem: $$\label{eq-scattering-quasiperiodic}
\left\{
\begin{array} {ll}
&{\displaystyle}\nabla \cdot \frac{1}{\rho_w} \nabla  v + \frac{\omega^2}{\kappa_w} v  = 0 \quad \text{in} \quad Y \backslash \overline{D}, \\
{\noalign{\smallskip}}&{\displaystyle}\nabla \cdot \frac{1}{\rho_b} \nabla  v+ \frac{\omega^2}{\kappa_b} v  = 0 \quad \text{in} \quad D, \\
{\noalign{\smallskip}}&{\displaystyle}v|_{+} -v|_{-}  =0   \quad \text{on} \quad \partial D, \\
{\noalign{\smallskip}}& {\displaystyle}\frac{1}{\rho_w} \frac{{\partial}v}{{\partial}\nu} \bigg|_{+} - \frac{1}{\rho_b} \frac{{\partial}v}{{\partial}\nu} \bigg|_{-} =0 \quad \text{on} \quad \partial D,\\
{\noalign{\smallskip}}&  e^{-i \alpha \cdot x} v  \,\,\,  \mbox{is periodic.}
  \end{array}
 \right.$$

It is known that has non-trivial solutions for discrete values of $\omega$ such as $$0 \le \omega_1^\alpha \le \omega_2^\alpha \le \cdots$$ and thus we have the following band structure of propagating frequencies for the periodic bubbly crystal ${\mathcal{C}}$: $$[0,\max_\alpha \omega_1^\alpha] \cup [\min_\alpha \omega_2^\alpha, \max_\alpha \omega_2^\alpha ] \cup \cdots.$$

In [@bandgap], 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 of which has the following form: $$\label{Helm-solution}
v =
\begin{cases}
\mathcal{S}_{D}^{\alpha,k_w} [\varphi^\alpha]\quad & \text{in} ~ Y \setminus \bar{D},\\
 \mathcal{S}_{D}^{k_b} [\psi^\alpha]   &\text{in} ~   {D},
\end{cases}$$ for some densities $\varphi^\alpha, \psi^\alpha \in  L^2({\partial}D)$. Using the jump relations for the single layer potentials, one can show that  is equivalent to the boundary integral equation $$\label{eq-boundary}
\mathcal{A}^\alpha(\omega, \delta)[\Phi^\alpha] =0,$$ where $$\mathcal{A}^\alpha(\omega, \delta) = 
 \begin{pmatrix}
  \mathcal{S}_D^{k_b} &  -\mathcal{S}_D^{\alpha,k}  \\
  -\frac{1}{2}+ \mathcal{K}_D^{k_b, *}& -\delta( \frac{1}{2}+ (\mathcal{K}_D^{ -\alpha,k})^*)
\end{pmatrix}, 
\,\, \Phi^\alpha= 
\begin{pmatrix}
\varphi^\alpha\\
\psi^\alpha
\end{pmatrix}.$$ Since it can be shown that $\omega=0$ is a characteristic value for the operator-valued analytic function $\mathcal{A}(\omega,0)$, we can conclude the following result by the Gohberg-Sigal theory [@lecturenotes; @Gohberg1971].

For any $\delta$ sufficiently small, there exists a characteristic value $\omega_1^\alpha= \omega_1^\alpha(\delta)$ to the operator-valued analytic function $\mathcal{A}^\alpha(\omega, \delta)$ such that $\omega_1^\alpha(0)=0$ and $\omega_1^\alpha$ depends on $\delta$ continuously.

The next theorem gives the asymptotic expansion of $\omega_1^\alpha$ as $\delta\rightarrow 0$.

[[[@bandgap]]{}]{} \[approx\_thm\] For $\alpha \ne 0$ and sufficiently small $\delta$, we have $$\begin{aligned}
\omega_1^\alpha= \sqrt{\frac{\delta \capacity_{D,\alpha}}{\pi R^2}} + O(\delta^{3/2}), \label{o_1_alpha}\end{aligned}$$ where the constant $\capacity_{D,\alpha}$ is given by $$\capacity_{D,\alpha}:= - \langle(\mathcal{S}_D^{\alpha,0})^{-1} [\chi_{\partial D}], \chi_{\partial D}\rangle.$$

Let $\omega_1^*=\max_\alpha \omega_1^\alpha$. The following theorem specifies the subwavelength bandgap opening.

[[[@bandgap]]{}]{} \[main\_bandgap\] For every ${\varepsilon}>0$, there exists $\delta_0>0$ and $\tilde \omega > \omega_1^*$ such that $$[ \omega_1^*+{\varepsilon}, \tilde\omega ] \subset [\max_\alpha \omega_1^\alpha, \min_\alpha \omega_2^\alpha]$$ for $\delta<\delta_0$.

Integral representation for bubbly crystals with a defect {#sec-2}
=========================================================

Bubbly crystals with a defect: problem formulation
--------------------------------------------------

(0,0) circle (10pt) node[$D_d$]{}; (0,0) circle (6pt) node\[yshift=-15pt, xshift=13pt \][$D$]{}; (1,0) circle (10pt); (0,1) circle (10pt); (1,1) circle (10pt); (-1,0) circle (10pt); (0,-1) circle (10pt); (1,-1) circle (10pt); (-1,1) circle (10pt); (-1,-1) circle (10pt); (1.6,0) node[$\cdots$]{}; (-1.55,0) node[$\cdots$]{}; (0,1.6) node[$\vdots$]{}; (0,-1.5) node[$\vdots$]{}; (1.7,1) node\[xshift=2pt\][$\rho_w, \kappa_w$]{}; (1,1) node[$\rho_b, \kappa_b$]{}; (0.5,0.5) – (0.5,-0.5) node\[yshift=-2pt, xshift=5pt \][$Y$]{} – (-0.5,-0.5) – (-0.5,0.5) – cycle;

Consider now a perturbed crystal, where the central disk $D$ is replaced by a defect disk $D_d$ of radius $R_d$ with $R_d<1/2$. We will consider both the cases when $R_d > R$ and $R_d < R$. Let ${\mathcal{C}}_d = D_d \cup \left( \cup_{n\in{\mathbb{Z}}^2\setminus\{0,0\}} D+n \right)$ be the perturbed crystal and let ${\varepsilon}= R_d-R \in (-R, 0) \cup (0,1/2-R)$ be the perturbation of the radius. We consider the following problem: $$\label{eq:scattering}
\left\{
\begin{array} {ll}
    &{\displaystyle}\nabla \cdot \frac{1}{\rho_w} \nabla  u+ \frac{\omega^2}{\kappa_w} u  = 0 \quad \text{in} \quad {\mathbb{R}}^2 \backslash {\mathcal{C}}_d, \\
    {\noalign{\smallskip}}&{\displaystyle}\nabla \cdot \frac{1}{\rho_b} \nabla  u+ \frac{\omega^2}{\kappa_b} u  = 0 \quad \text{in} \quad {\mathcal{C}}_d, \\
    {\noalign{\smallskip}}&{\displaystyle}u|_{+} -u|_{-}  =0   \quad \text{on} \quad \partial {\mathcal{C}}_d, \\
    {\noalign{\smallskip}}& {\displaystyle}\frac{1}{\rho_w} \frac{\partial u}{\partial \nu} \bigg|_{+} - \frac{1}{\rho_b} \frac{\partial u}{\partial \nu} \bigg|_{-} =0 \quad \text{on} \quad \partial {\mathcal{C}}_d .
\end{array}
\right.$$ As discussed in Subsection \[subsec:bandgap\], the unperturbed problem (with ${\mathcal{C}}$ instead of ${\mathcal{C}}_d$ in [(\[eq:scattering\])]{}), has the following essential spectrum for the propagating frequencies: $$[0,\max_\alpha \omega_1^\alpha] \cup [\min_\alpha \omega_2^\alpha, \max_\alpha \omega_2^\alpha ] \cup \cdots.$$ In fact, it can be easily shown that the perturbed crystal problem [(\[eq:scattering\])]{} has the same essential spectrum. This is because the essential spectrum is stable under compact perturbations [@figotin; @MMMP4].

In this paper, we want to show that by modifying the central bubble $D$, there exists a frequency $\omega^{\varepsilon}$, slightly above $\max_\alpha\omega_1^\alpha$, which results in a non-trivial solution to the problem [(\[eq:scattering\])]{}. The solution $u$ associated with the frequency $\omega^{\varepsilon}$ should be localized since $\omega^{\varepsilon}$ lies in the bandgap. Moreover, it is of a subwavelength nature.

Effective sources for the defect
--------------------------------

It is difficult to obtain the boundary integral formulation of the problem directly. Here we consider the effective source solution which models the defect $D_d$ by placing non-trivial sources onto the boundary of the central bubble $D$ of the unperturbed crystal ${\mathcal{C}}$. Since the geometry of the unperturbed crystal is periodic, we can use the Floquet-Bloch theory. This is motivated by the fictitious source superposition method introduced in [@Wilcox].

Let us consider the following problem: $$\label{eq:scattering_fictitious}
\left\{
\begin{array} {ll}
    &{\displaystyle}\nabla \cdot \frac{1}{\rho_w} \nabla  \tilde{u}+ \frac{\omega^2}{\kappa_w} \tilde{u}  = 0 \quad \text{in} \quad {\mathbb{R}}^2 \setminus {\mathcal{C}}, \\
    {\noalign{\smallskip}}&{\displaystyle}\nabla \cdot \frac{1}{\rho_b} \nabla  \tilde{u}+ \frac{\omega^2}{\kappa_b} \tilde{u}  = 0 \quad \text{in} \quad {\mathcal{C}}, \\
    {\noalign{\smallskip}}&{\displaystyle}\tilde{u}|_{+} -\tilde{u}|_{-}  =f\delta_{m,(0,0)}  \quad \text{on} \quad \partial D+m, \ m\in\mathbb{Z}^2, \\
    {\noalign{\smallskip}}& {\displaystyle}\frac{1}{\rho_w} \frac{\partial \tilde{u}}{\partial \nu} \bigg|_{+} - \frac{1}{\rho_b} \frac{\partial \tilde{u}}{\partial \nu} \bigg|_{-} =g\delta_{m,(0,0)} \quad \text{on} \quad \partial D+m, \ m\in\mathbb{Z}^2.
\end{array}
\right.$$ where $f,g$ are the source terms and $\delta_{m,n}$ is the Kronecker delta function. Note that the non-zero sources are present only on the boundary of the central bubble $D$.

Introduce the notation $$D_l = \begin{cases} D_d \quad & \text{if } {\varepsilon}> 0 \\ D \quad & \text{if } {\varepsilon}< 0, \end{cases} \qquad
D_s = \begin{cases} D \quad & \text{if } {\varepsilon}> 0 \\ D_d \quad & \text{if } {\varepsilon}< 0. \end{cases}$$ Thus, $D_l$ and $D_s$ denotes the largest and smallest of $D, D_d$, respectively. If the source terms $f$ and $g$ satisfy some appropriate conditions, then we will have $$u \equiv \tilde{u}  \quad \mbox{in } (\mathbb{R}^2\setminus D_l) \cup D_s.$$ Once this is achieved, the original solution $u$ can be recovered by extending the solution $\widetilde{u}$ to the whole region including $D_l\setminus D_s$ with boundary conditions on ${\partial}D$ and ${\partial}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 {#subsec:char_eff}
-----------------------------------------

Here we clarify the relation between the effective source pair $(f,g)$ and the density pair $(\varphi,\psi)$.

First we consider the integral equation for the solution $\widetilde{u}$ inside the central unit cell $Y$. Inside $Y$, the solution $\tilde{u}$ can be decomposed as $$\begin{aligned}
\label{eq:u_rep_effective}
\tilde{u}
= 
\begin{cases}
H + {\mathcal{S}}_D^{k_w}[\psi] &\quad \mbox{in } Y\setminus \overline{D},
\\
{\mathcal{S}}_D^{k_b}[\varphi] &\quad \mbox{in } D,
\end{cases}\end{aligned}$$ where $H$ satisfies the homogeneous equation $(\Delta + k_w^2) H = 0$ in $Y$ and the pair $(\varphi,\psi)\in L^2({\partial}D)^2$ satisfies $$\label{eq:AD}
\mathcal{A}_D
\begin{pmatrix}
\varphi 
\\
\psi
\end{pmatrix}
:=
\begin{pmatrix}
{\mathcal{S}}_D^{k_b} & - {\mathcal{S}}_D^{k_w}
\\[0.3em]
{\displaystyle}\partial_\nu {\mathcal{S}}_D^{k_b} |_{-}
&
{\displaystyle}-\delta \partial_\nu {\mathcal{S}}_D^{k_w} |_{+}
\end{pmatrix}
\begin{pmatrix}
\varphi 
\\
\psi
\end{pmatrix}
=
\begin{pmatrix}
H|_{\partial D}- f
\\[0.3em]
{\displaystyle}\partial_\nu H |_{\partial D} - g
\end{pmatrix}.$$

Next, let us consider the central cell $Y$ in the original problem . The central cell $Y$ contains a defect bubble $D_d$ and no sources. Inside $Y$, the solution $u$ is represented as $$\begin{aligned}
u=
\begin{cases}
H + {\mathcal{S}}_{D_d}^{k_w}[\psi_d] &\quad \mbox{in } Y\setminus \overline{D_d},
\\
{\mathcal{S}}_{D_d}^{k_b}[\varphi_d] &\quad \mbox{in } D_d,
\end{cases}\end{aligned}$$ where $$\label{eq:ADd_original}
{\mathcal{A}}_{D_d}
\begin{pmatrix}
\varphi_d 
\\
\psi_d
\end{pmatrix}:=
\begin{pmatrix}
{\mathcal{S}}_{D_d}^{k_b} & - {\mathcal{S}}_{D_d}^{k_w}
\\[0.3em]
{\displaystyle}\partial_\nu {\mathcal{S}}_{D_d}^{k_b} |_{-}
&
{\displaystyle}-\delta \partial_\nu {\mathcal{S}}_{D_d}^{k_w} |_{+}
\end{pmatrix}
\begin{pmatrix}
\varphi_d 
\\
\psi_d
\end{pmatrix}
=
\begin{pmatrix}
H|_{\partial {D_d}}
\\[0.3em]
{\displaystyle}\partial_\nu H |_{\partial {D_d}}
\end{pmatrix}.$$

In order for the effective sources $f$ and $g$ to model the defect $D_d$ correctly, we should impose $u \equiv \widetilde{u} \   \mbox{in } (Y\setminus D_l) \cup D_s$, which implies $u \equiv \widetilde{u} \  \mbox{in } \mathbb{R}^2\setminus Y $ as well by an argument using analytic continuation. In other words, the following conditions should be satisfied: $$\label{eq:SDdSD_inside}
\mathcal{S}_{D_d}^{k_b}[\varphi_d] \equiv \mathcal{S}_D^{k_b}[\varphi] \quad\mbox{in }D_s,$$ and $$\label{eq:SDdSD_outside}
\mathcal{S}_{D_d}^{k_w}[\psi_d] \equiv \mathcal{S}_D^{k_w}[\psi] \quad\mbox{in }Y\setminus \overline{D_l}.$$

Since $D$ and $D_d$ are circular disks, we can use a Fourier basis for functions on ${\partial}D$ or ${\partial}D_d$ in polar coordinates $(r,\theta)$ to make - more explicit. Let us write $\varphi$ and $\psi$ in the form of Fourier series: $$\varphi(\theta) = \sum_{n\in\mathbb{Z}} \varphi_n e^{in \theta}, \quad
\psi(\theta) = \sum_{n\in\mathbb{Z}} \psi_n e^{in \theta}.$$ Similarly, we also write $\varphi_d$ and $\psi_d$ as $$\varphi_d(\theta) = \sum_{n\in\mathbb{Z}} \varphi_{d,n} e^{in \theta}, \quad
\psi_d(\theta) = \sum_{n\in\mathbb{Z}} \psi_{d,n} e^{in \theta}.$$ We define the subspace $V_{mn}$ of $L^2({\partial}D)^2$ as $$V_{mn} : = \mbox{span}\{ e^{im \theta}\} \times \mbox{span}\{ e^{in \theta}\}, \quad
m,n\in \mathbb{Z}.$$ Similarly, let ${V}_{mn}^d$ be the subspace of $L^2({\partial}D_d)^2$ with the same Fourier basis.

It is known that (see, for instance, [@bandgap]) $$\label{eq:single_explicit}
{\mathcal{S}}_D^k[e^{i n \theta}]
=\frac{(-1)i\pi R}{2}\times \begin{cases}
J_n(k R) H_n^{(1)}(k r) e^{i n\theta}, &\quad r\geq R,
\\
H_n^{(1)}(k R) J_n(k r) e^{i n\theta}, &\quad 0\leq r < R. 
\end{cases}$$ Therefore, the operator ${\mathcal{A}}_{D}$ in has the following matrix representation as an operator from $V_{mn}$ to $V_{m'n'}$: $$\label{eq:AD_multipole}
({\mathcal{A}}_{D})_{V_{mn}\rightarrow V_{m'n'}} = \delta_{mn}\delta_{m'n'} \frac{(-i)\pi R}{2} \begin{pmatrix}
  J_n(k_b R)H_n^{(1)}(k_b R) & -  J_n(k_w R)H_n^{(1)}(k_w R)
 \\
 k_b J_n'(k_b R) H_n^{(1)}(k_b R) & - \delta k_w J_n(k_w R) \big(H_n^{(1)}\big)'(k_w R)
 \end{pmatrix}.$$ Similarly, the operator ${\mathcal{A}}_{D_d}$ in , as a mapping from $V_{mn}^d$ to $V_{m'n'}^d$, is represented as follows: $$\label{eq:ADd_multipole}
({\mathcal{A}}_{D_d})_{V_{mn}^d\rightarrow V_{m'n'}^d} = \delta_{mn}\delta_{m'n'} \frac{(-i)\pi R_d}{2} \begin{pmatrix}
  J_n(k_b R_d)H_n^{(1)}(k_b R_d) & -  J_n(k_w R_d)H_n^{(1)}(k_w R_d)
 \\
 k_b J_n'(k_b R_d) H_n^{(1)}(k_b R_d) & - \delta k_w J_n(k_w R_d) \big(H_n^{(1)}\big)'(k_w R_d)
 \end{pmatrix}.$$

Now, we consider and . We have from that, inside $D_s$, $$\begin{aligned}
\mathcal{S}_{D}^{k_b}[\varphi] &= \sum_{n\in\mathbb{Z}}\varphi_{n} \mathcal{S}_{D}^{k_b}[e^{i n \theta}]
= \frac{(-i)\pi R}{2}\sum_{n\in\mathbb{Z}}\varphi_{n}  H_n^{(1)}(k_b R) J_n(k_b r) e^{i n\theta}
\\
&= \frac{(-i)\pi R}{2}\sum_{n\in\mathbb{Z}}\varphi_{n}  \frac{H_n^{(1)}(k_b R)}{H_n^{(1)}(k_b R_d)}  H_n^{(1)}(k_b R) J_n(k_b r) e^{i n\theta}
\\
&= \mathcal{S}_{D_d}^{k_b}\Big[\sum_{n\in\mathbb{Z}}\varphi_{n} \frac{R}{R_d}\frac{H_n^{(1)}(k_b R)}{H_n^{(1)}(k_b R_d)}e^{in\theta}\Big].\end{aligned}$$ Similarly, outside $D_l$, $$\mathcal{S}_{D}^{k_w}[\psi] = \mathcal{S}_{D_d}^{k_w}\Big[\sum_{n\in\mathbb{Z}}\psi_{n} \frac{R}{R_d} \frac{J_n(k_w R)}{J_n(k_w R_d)}e^{in\theta}\Big].$$ So, from and , we see that $$\begin{pmatrix}
\varphi_d
\\[0.3em]
{\displaystyle}\psi_d
\end{pmatrix}
=
\mathcal{P}_1
\begin{pmatrix}
\varphi
\\[0.3em]
\psi
\end{pmatrix},$$ where the operator $\mathcal{P}_1: L^2({\partial}D)^2\rightarrow L^2({\partial}D_d)^2$ is given by $$(\mathcal{P}_1)_{V_{mn}\rightarrow V_{m'n'}^d} = \delta_{mn}   \delta_{m'n'}\frac{R}{R_d}
\begin{pmatrix}
{\displaystyle}\frac{H_n^{(1)}(k_b R) }{H_n^{(1)}(k_b R_d)}
& 0
\\
0& {\displaystyle}\frac{J_n(k_w R) }{J_n(k_w R_d)}
\end{pmatrix}.
\quad$$ In the same way, we can obtain that $$\begin{pmatrix}
H|_{\partial D_d}
\\[0.3em]
{\displaystyle}\partial_\nu H |_{\partial D_d}
\end{pmatrix}
=
\mathcal{P}_2
\begin{pmatrix}
H|_{\partial D}
\\[0.3em]
{\displaystyle}\partial_\nu H |_{\partial D}
\end{pmatrix},$$ where the operator $\mathcal{P}_2: L^2({\partial}D)^2\rightarrow L^2({\partial}D_d)^2 $ is given by $$(\mathcal{P}_2)_{V_n\rightarrow V_m^d} = \delta_{mn}  
\begin{pmatrix}
{\displaystyle}\frac{J_n(k_w R_d) }{J_n(k_w R)}
& 0
\\
0& {\displaystyle}\frac{J_n'(k_w R_d) }{J_n'(k_w R)}
\end{pmatrix}.$$ Therefore, using , we arrive at $$\label{eq:ADd}
\mathcal{A}_D^{\varepsilon}\begin{pmatrix}
\varphi
\\[0.3em]
\psi
\end{pmatrix} := (\mathcal{P}_2)^{-1}\mathcal{A}_{D_d} \mathcal{P}_1 
\begin{pmatrix}
\varphi
\\[0.3em]
\psi
\end{pmatrix}
=\begin{pmatrix}
H|_{\partial D}
\\[0.3em]
{\displaystyle}\partial_\nu H |_{\partial D}
\end{pmatrix}.$$ Thus, yields $$\label{eq:relation_density_source}
(\mathcal{A}_D^{\varepsilon}- \mathcal{A}_D)\begin{pmatrix}
\varphi
\\[0.3em]
\psi
\end{pmatrix} = 
\begin{pmatrix}
f
\\[0.3em]
g
\end{pmatrix}.$$

We have obtained an explicit relation between the pair $(\varphi,\psi)$ and the effective sources $(f,g)$. If the effective sources $f$ and $g$ satisfy , then the pair $(f,g)$ will result in the generation of the same scattered field as the one induced by the defect bubble $D_d$. In other words, $\widetilde{u}\equiv u$ outside $D_l\setminus \overline{D_s}$. In what follows, for convenience of notation, we will identify the solution $\widetilde u$ with the original one $u$.

Floquet transform of the solution
---------------------------------

Here we derive an integral equation for the effective source problem . We apply the Floquet transform to the solution $u$ with the quasi-periodic parameter $\alpha$ as follows $$u^\alpha = \sum_{m\in {\mathbb{Z}}^2} u(x-m) e^{i\alpha \cdot m}.$$ The transformed solution $u^\alpha$ satisfies $$\label{eq:scattering_quasi}
\left\{
\begin{array} {ll}
    &{\displaystyle}\nabla \cdot \frac{1}{\rho_w} \nabla  u^\alpha+ \frac{\omega^2}{\kappa_w} u^\alpha  = 0 \quad \text{in} \quad Y \setminus \overline{D}, \\
    {\noalign{\smallskip}}&{\displaystyle}\nabla \cdot \frac{1}{\rho_b} \nabla  u^\alpha+ \frac{\omega^2}{\kappa_b} u^\alpha  = 0 \quad \text{in} \quad D, \\
    {\noalign{\smallskip}}&{\displaystyle}u^\alpha|_{+} -u^\alpha|_{-}  =f   \quad \text{on} \quad \partial D, \\
    {\noalign{\smallskip}}& {\displaystyle}\frac{1}{\rho_w} \frac{\partial u^\alpha}{\partial \nu} \bigg|_{+} - \frac{1}{\rho_b} \frac{\partial u^\alpha}{\partial \nu} \bigg|_{-} =g \quad \text{on} \quad \partial D,
    \\
    {\noalign{\smallskip}}& {\displaystyle}e^{-i \alpha \cdot x} u^\alpha \text{ is periodic}.
\end{array}
\right.$$ The solution $u^\alpha$ can be represented using quasi-periodic layer potentials as $$\begin{aligned}
u^\alpha
= 
\begin{cases}
{\mathcal{S}}_D^{\alpha,k_w}[\psi_\alpha], &\quad \mbox{in } Y\setminus \overline{D},
\\
{\mathcal{S}}_D^{k_b}[\varphi_\alpha], &\quad \mbox{in } D,
\end{cases}\end{aligned}$$ where the pair $(\varphi^\alpha,\psi^\alpha)\in L^2({\partial}D)^2$ is the solution to $$\label{eq:phipsialpha}
\mathcal{A}^\alpha(\omega,\delta)
\begin{pmatrix}
\varphi^\alpha 
\\[0.3em]
\psi^\alpha
\end{pmatrix}
:=
\begin{pmatrix}
{\mathcal{S}}_D^{k_b} & - {\mathcal{S}}_D^{\alpha,k_w}
\\[0.3em]
{\displaystyle}\partial_\nu {\mathcal{S}}_D^{k_b} |_{-}
&
{\displaystyle}-\delta \partial_\nu {\mathcal{S}}_D^{\alpha,k_w} |_{+}
\end{pmatrix}
\begin{pmatrix}
\varphi^\alpha 
\\[0.3em]
\psi^\alpha
\end{pmatrix}
=
\begin{pmatrix}
- f
\\[0.3em]
- g
\end{pmatrix}.$$ Since the operator $\mathcal{A}_\alpha$ is invertible for $\omega$ in the bandgap [@bandgap], we have $$\begin{pmatrix}
\varphi^\alpha 
\\[0.3em]
\psi^\alpha
\end{pmatrix} 
= \mathcal{A}^\alpha(\omega,\delta)^{-1}\begin{pmatrix}
- f
\\[0.3em]
- g
\end{pmatrix}.$$

The original solution $u$ can be recovered by the inversion formula as follows $$u(x)=\frac{1}{(2\pi)^2}\int_{BZ} u^\alpha(x) d\alpha.$$ Then, inside the region $D$, the solution $u$ satisfies $$u= {\mathcal{S}}_D^{k_b}\big[\frac{1}{(2\pi)^2}\int_{BZ}\varphi^\alpha d\alpha \big].$$ Similarly, inside the region $Y\setminus \overline{D}$, we have $$\begin{aligned}
{\displaystyle}u&=\frac{1}{(2\pi)^2}\int_{BZ}{\mathcal{S}}_D^{\alpha,k_w}[\psi^\alpha] d\alpha
={\mathcal{S}}_D^{k_w}\big[\frac{1}{(2\pi)^2}\int_{BZ}\psi^\alpha  d\alpha\big] + \frac{1}{(2\pi)^2}\int_{BZ} \sum_{m\in\mathbb{Z}^2,m\neq 0}{\mathcal{S}}_D^{k_w}[\psi^\alpha](\,\cdot - m)e^{im\cdot \alpha} d\alpha.\end{aligned}$$ Note that the second term in the right-hand side satisfies the homogeneous Helmholtz equation $(\Delta + k_w^2 )u=0$ in $Y\setminus \overline{D}$. So, in view of , we can identify $\varphi,\psi$ and $H$ as follows: $$\varphi = \frac{1}{(2\pi)^2}\int_{BZ} \varphi^\alpha d\alpha,
\quad
\psi = \frac{1}{(2\pi)^2}\int_{BZ} \psi^\alpha d \alpha,$$ and $$H=\frac{1}{(2\pi)^2}\int_{BZ} \sum_{m\in\mathbb{Z}^2,m\neq 0}{\mathcal{S}}_D^{k_w}[\psi^\alpha](\,\cdot - m)e^{im\cdot \alpha} d\alpha.$$ Therefore, from , we get the following result.

The density pair $(\varphi,\psi)$ and the effective source pair $(f,g)$ satisfy $$\label{eq:IntAa}
\begin{pmatrix}
\varphi
\\[0.3em]
\psi 
\end{pmatrix}
=\left(\frac{1}{(2\pi)^2}\int_{BZ} \mathcal{A}^\alpha(\omega,\delta)^{-1} d\alpha \right)\begin{pmatrix}
-f
\\[0.3em]
-g 
\end{pmatrix}$$ for $\omega$ in the bandgap.

The integral equation for the effective sources
-----------------------------------------------

Here we derive the integral equation for the effective source pair $(f,g)$. We have the following result.

The effective source pair $(f,g)\in L^2(\partial D)^2$ satisfies the following integral equation: $$\begin{aligned}
\label{eq:fg}
\mathcal{M}^{\varepsilon}(\omega)\begin{pmatrix}
f
\\
g
\end{pmatrix}:=\bigg(I+({\mathcal{A}}_{D}^{\varepsilon}(\omega,\delta) - {\mathcal{A}}_{D}(\omega,\delta))\frac{1}{(2\pi)^2}\int_{BZ}{\mathcal{A}}^\alpha(\omega,\delta)^{-1} d\alpha\bigg)
\begin{pmatrix}
f 
\\[0.3em]
 g
\end{pmatrix}
=
\begin{pmatrix}
0
\\[0.3em]
0
\end{pmatrix},\end{aligned}$$ for $\omega$ in the bandgap.

Equation [(\[eq:fg\])]{} is an immediate consequence of equations [(\[eq:IntAa\])]{} and [(\[eq:relation\_density\_source\])]{}.

Thus, if we find a value of $\omega$ in the bandgap such that there exists a non-trivial solution pair $(f,g)$ to (\[eq:fg\]), then we will find a resonant frequency for the localized mode.

Subwavelength localized modes {#sec-3}
=============================

The resonant frequency of the localized mode
--------------------------------------------

Here we prove that a frequency for the localized mode exists slightly above $\max_\alpha \omega_1^\alpha$. In what follows, let us omit the subscript in $\omega_1^\alpha$ for ease of notation. We also do not make explicit the dependence on $\delta$.

We need to study the characteristic value of the operator $\mathcal{M}^{\varepsilon}$ appearing in . Let us first analyse the operator $\int_{BZ} ({\mathcal{A}}^{\alpha})^{-1}d\alpha$. Since $\omega^\alpha$ is a simple pole of the mapping $\omega \mapsto {\mathcal{A}}^\alpha(\omega,)^{-1}$ in a neighbourhood of $\omega^\alpha$, according to [@lecturenotes], we can write $${\mathcal{A}}^\alpha(\omega)^{-1} = \frac{{\mathcal{L}}^\alpha}{\omega- \omega^\alpha} + \mathcal{R}^\alpha(\omega),$$ where the operator-valued function $\mathcal{R}^\alpha(\omega)$ is holomorphic in a neighbourhood of $\omega^\alpha$, and the operator ${\mathcal{L}}^\alpha$ maps $L^2({\partial}D)^2$ onto $\ker {\mathcal{A}}^\alpha(\omega^\alpha,\delta)$. Let us write $$\ker {\mathcal{A}^\alpha(\omega^\alpha)} = \mbox{span} \{\Psi^\alpha\}, \quad \ker {\big(\mathcal{A}^\alpha(\omega^\alpha)\big)^\dagger} = \mbox{span} \{\Phi^\alpha\},$$ where $^\dagger$ denotes the adjoint operator. Then, as in [@lecturenotes; @thinlayer], it can be shown that $${\mathcal{L}}^\alpha = \frac{\langle \Phi^\alpha,\ \cdot\ \rangle \Psi^\alpha}{\langle \Phi^\alpha, \frac{d}{d \omega}{\mathcal{A}}^\alpha\big|_{\omega=\omega^\alpha} \Psi^\alpha \rangle},$$ where $\langle \,\cdot \,,\,\cdot\, \rangle$ stands for the standard inner product of $L^2(\partial D)^2$.

Hence the operator $\mathcal{M}^{\varepsilon}$ can be decomposed as $$\mathcal{M}^{\varepsilon}= I + ({\mathcal{A}}_{D}^{\varepsilon}- {\mathcal{A}}_D) \frac{1}{(2\pi)^2}\int_{BZ}\frac{{\mathcal{L}}^\alpha}{\omega-\omega_\alpha}d\alpha +  ({\mathcal{A}}_{D}^{\varepsilon}- {\mathcal{A}}_D) \frac{1}{(2\pi)^2}\int_{BZ} \mathcal{R}_\alpha d\alpha.$$ Note that the third term in the right-hand side is holomorphic with respect to $\omega$.

Denote by $\alpha^* = (\pi,\pi)$ and $\omega^* = \omega^{(\pi,\pi)}$. In [@highfrequency], it was proved that, using the symmetry of the square array of bubbles, $\omega^\alpha$ attains its maximum at $\alpha=\alpha^*$. Since we are assuming that each bubble is a circular disk, we can derive a slightly more refined result as shown in the following lemma.

\[lem:max\_omegaalp\] The characteristic value $\omega_\alpha$ attains its maximum at $\alpha=\alpha^*$. Moreover, for $\alpha$ near $\alpha_*$, we have $$\omega^\alpha = \omega^* -   \frac{1}{2}c_\delta |\alpha-\alpha_*|^2   + o(|\alpha-\alpha_*|^2).$$ Here, $c_\delta$ is a positive constant depending on $\delta$ and scales as $c_\delta = O(\sqrt{\delta})$.

In view of the above lemma, we can expect that the operator $\int_{BZ}\frac{{\mathcal{L}}^\alpha}{\omega-\omega^\alpha}$ becomes singular when $\omega\rightarrow \omega^*(=\max_{\alpha } \omega^{\alpha})$. Let us extract its singular part explicitly. Before doing this, we introduce some notations. Denote by $
{\mathcal{A}}^* = {\mathcal{A}}^{(\pi,\pi)},  \Phi^* = \Phi^{(\pi,\pi)},$ and ${\mathcal{L}}^* = {\mathcal{L}}^{(\pi,\pi)}
$. We also define a small neighbourhood $V$ of $\omega^*$ which excludes the real interval $(-\infty,\omega^*]$, namely, $$V= \{|\omega-\omega^*|<r_*\} \setminus (-\infty,\omega^*],$$ for some small enough $r_*>0$. By Lemma \[lem:max\_omegaalp\], we have $$\begin{aligned}
\frac{1}{(2\pi)^2}\int_{BZ}\frac{{\mathcal{L}}^\alpha}{\omega-\omega^\alpha} d\alpha &= {\mathcal{L}}^*\frac{1}{(2\pi)^2}\int_{BZ} \frac{1}{\omega - \omega^*+\frac{1}{2} c_\delta |\alpha-\alpha^*|^2} d\alpha + B_1(\omega)
\\
&={\mathcal{L}}^*\frac{1}{(2\pi)^2}\int_{|\alpha-\alpha^*|<1} \frac{}{\omega - \omega^*+\frac{1}{2} c_\delta |\alpha-\alpha^*|^2} d\alpha + B_2(\omega).\end{aligned}$$ Hereafter, $B_j$ means a bounded function with respect to $\omega$ in $V$. Then, using the polar coordinates $\alpha-\alpha_*=(r',\theta')$, we get $$\begin{aligned}
\frac{1}{(2\pi)^2}\int_{BZ}\frac{{\mathcal{L}}^\alpha}{\omega-\omega^\alpha} d\alpha
&=\frac{1}{(2\pi)^2}\int_0^{1} \frac{2\pi r'}{\omega-\omega_* + \frac{1}{2} c_\delta r'^2} dr' +B_2(\omega)
\\
&= -\frac{1}{2\pi c_\delta}\ln (\omega-\omega_*) + B_3(\omega).\end{aligned}$$ Here, the usual principal branch is taken for the logarithm and so the operator $\int_{BZ}\frac{{\mathcal{L}}^\alpha}{\omega-\omega^\alpha} d\alpha$ has a branch cut on $(-\infty,\omega^*]$. We also observe that, for $\omega \in V$, $$\ln(\omega-\omega^*) ({\mathcal{A}}_{D}^{\varepsilon}(\omega)-{\mathcal{A}}_{D}(\omega))=
\ln(\omega-\omega^*) ({\mathcal{A}}_{D}^{\varepsilon}(\omega^*)-{\mathcal{A}}_{D}(\omega)^*) + B_4(\omega).$$

Therefore, the integral equation can be rewritten as $$\begin{aligned}
\mathcal{M}^{\varepsilon}(\omega)\begin{pmatrix}
f
\\
g
\end{pmatrix}= \Big(I - \frac{\ln(\omega-\omega^*)}{2\pi c_\delta} ({\mathcal{A}}_{D}^{\varepsilon}(\omega^*)-{\mathcal{A}}_{D}(\omega^*))\mathcal{L}^*  + {\mathcal{R}^{\varepsilon}(\omega)}
\Big)
\begin{pmatrix}
f
\\
g
\end{pmatrix}
=
\begin{pmatrix}
0
\\
0
\end{pmatrix},\end{aligned}$$ where $\mathcal{R}^{\varepsilon}$ is analytic and bounded with respect to $\omega$ in $V$ and satisfies $
{\mathcal{R}}^{\varepsilon}= O({\varepsilon})
$. We first consider the principal part $\mathcal{N}^{\varepsilon}(\omega)$, namely, $$\mathcal{N}^{\varepsilon}: \omega \mapsto \mathcal{N}^{\varepsilon}(\omega)=
I - \frac{1}{2\pi c_\delta}\ln(\omega-\omega^*) ({\mathcal{A}}_{D}^{\varepsilon}(\omega^*)-{\mathcal{A}}_{D}(\omega^*))\mathcal{L}^*.$$ Let us find its characteristic value $\widehat\omega$, [*i.e.*]{}, $\widehat\omega\in V$ such that there exists a non-trivial function $\widehat\Phi$ satisfying $\mathcal{N}^{\varepsilon}(\widehat\omega) \widehat\Phi =0$. Equivalently, we have $$\begin{aligned}
\widehat\Phi - \frac{1}{2\pi c_\delta}\ln(\widehat\omega-\omega^*) \frac{ ({\mathcal{A}}_{D}^{\varepsilon}-{\mathcal{A}}_{D})(\omega^*) \Psi^*}{\langle \Phi^*, \frac{d}{d \omega}{\mathcal{A}}^*\big|_{\omega=\omega^*} \Psi^* \rangle}\langle \Phi^*, \widehat\Phi \rangle =0.\end{aligned}$$ Then, by multiplying by $\Phi^*$, $$\begin{aligned}
\langle \Phi^*, \widehat\Phi\rangle \bigg( 1 -  \frac{1}{2\pi c_\delta}\ln(\widehat\omega-\omega^*) \frac{ \langle \Phi^*,({\mathcal{A}}_{D}^{\varepsilon}-{\mathcal{A}}_{D})(\omega^*) \Psi^*\rangle}{\langle \Phi^*, \frac{d}{d \omega}{\mathcal{A}}^*\big|_{\omega=\omega^*} \Psi^* \rangle} \bigg)=0.\end{aligned}$$ Since $\langle \Phi^*, \widehat\Phi \rangle=0$ would imply $\widehat\Phi=0$, we get $$\label{eq:temptemp}
1 -  \frac{1}{2\pi c_\delta}\ln(\widehat\omega-\omega^*) \frac{ \langle \Phi^*,({\mathcal{A}}_{D}^{\varepsilon}-{\mathcal{A}}_{D})(\omega^*) \Psi^*\rangle}{\langle \Phi^*, \frac{d}{d \omega}{\mathcal{A}}^*\big|_{\omega=\omega^*} \Psi^* \rangle} = 0.$$

Before solving the above equation for $\omega$, we need the following Lemma whose proof will be given in Subsection \[subsec:proofs\_estim\].

\[lem:estim1\] The following results hold:

-   When $\delta\rightarrow 0$, we have $$\big\langle \Phi^*, \frac{d}{d \omega}{\mathcal{A}}^*(\omega^*,\delta) \Psi^* \big\rangle = -2\pi{\omega^*\ln \omega^*} R^3 + O(\sqrt{\delta}),$$ which is positive for $\delta$ small enough.

-   For a fixed ${\varepsilon}$, when $\delta\rightarrow 0$ we have $$\begin{aligned}
            \left\langle\Phi^*,\big({\mathcal{A}}_{D}^{\varepsilon}(\omega^*,\delta)-{\mathcal{A}}_{D}(\omega^*,\delta)\big){ \Psi^*} \right\rangle &=\delta{\varepsilon}\ln\omega^*\left(R\|\psi_{\alpha^*}\|^2_{L^2({\partial}D)}  - 2\capacity_{D,\alpha^*} \right) + O({\varepsilon}\delta + {\varepsilon}^2\delta\ln\delta),\\
            & := {\varepsilon}\big(\delta\ln\omega^*S(R)\big) + O({\varepsilon}\delta + {\varepsilon}^2\delta\ln\delta),
            \end{aligned}$$ where $\psi_{\alpha^*} = ({\mathcal{S}}_{D}^{\alpha^*,0})^{-1}[\chi_{{\partial}D}]$ and $\capacity_{D,\alpha^*} = -\langle \psi_{\alpha^*}, \chi_{{\partial}D}\rangle$. For small ${\varepsilon}$ and $\delta$, the sign of $\delta\ln\omega^*S(R)$ varies with $R$ as follows: $\delta\ln\omega^*S(R) > 0$ for small enough $R$, while $\delta\ln\omega^*S(R) < 0$ for $R$ close enough to $1/2$.

In view of Lemma \[lem:estim1\], we have two different regimes: The dilute regime when $R$ is small and the non-dilute regime when $R$ is close to $1/2$. For equation to have a real solution $\widehat\omega$ close to $\omega_*$, we need ${\varepsilon}< 0$ in the dilute regime and ${\varepsilon}> 0$ in the non-dilute regime. Under this assumption we have $$\begin{aligned}
\widehat\omega-\omega_* = 
  \exp\big(2\pi{c_\delta} \big\langle \Phi^*, \frac{d}{d \omega}{\mathcal{A}}^*(\omega^*) \Psi^* \big\rangle \big\langle\Phi^*,\big({\mathcal{A}}_{D}^{\varepsilon}(\omega^*)-{\mathcal{A}}_{D}(\omega^*)\big){ \Psi^*} \big\rangle^{-1}\big).\end{aligned}$$ We can also see that the right-hand side is positive and goes to zero as ${\varepsilon}$ tends to zero. In other words, $\widehat\omega\rightarrow \omega^*$ as ${\varepsilon}\rightarrow 0$. Now we turn to the full operator $\mathcal{M}^{\varepsilon}$. Recall that $\mathcal{M}^{\varepsilon}= \mathcal{N}^{\varepsilon}+ \mathcal{R}^{\varepsilon}$and $\mathcal{R}^{\varepsilon}$ is holomorphic, bounded and satisfies $\mathcal{R}^{\varepsilon}=O({\varepsilon})$ in $V$. So by Gohberg-Sigal theory [@lecturenotes], we can conclude that there exists a characteristic value $\omega^{\varepsilon}$ of the operator-valued function $\mathcal{M}^{\varepsilon}$ near $\omega^*$. Let us denote its associated root function by $\Phi^{\varepsilon}$. We choose $\Phi^{\varepsilon}$ so that $\langle \Phi^*, \Phi^{\varepsilon}\rangle=1$. Then, as in the derivations of , we can obtain $$1 - \frac{1}{2\pi c_\delta}\ln(\omega^{\varepsilon}-\omega^*) \mu({\varepsilon},\delta,R)  + \langle \Phi^*,\mathcal{R}^{\varepsilon}(\omega^{\varepsilon}) \Phi^{\varepsilon}\rangle= 0,$$ where $$\mu({\varepsilon},\delta,R) = \frac{ \langle \Phi^*,({\mathcal{A}}_{D}^{\varepsilon}-{\mathcal{A}}_{D})(\omega^*) \Psi^*\rangle}{\langle \Phi^*, \frac{d}{d \omega}{\mathcal{A}}^*\big|_{\omega=\omega^*} \Psi^* \rangle}.$$ Note that $\langle \Phi^*,\mathcal{R}^{\varepsilon}(\omega) \Phi^{\varepsilon}\rangle = O({\varepsilon})$ for small ${\varepsilon}$. We see from Lemma \[lem:estim1\] that $${\mu}({\varepsilon},\delta,R) =  \frac{\delta{\varepsilon}\left(R\|\psi_{\alpha^*}\|^2_{L^2({\partial}D)}  - 2\capacity_{D,\alpha^*} \right)}{-2\pi\omega^* R^3} + O\Big({\varepsilon}\frac{\sqrt{\delta}}{\ln\delta} + \sqrt{\delta}{\varepsilon}^2 \Big),$$ as $\delta\rightarrow 0$ and ${\varepsilon}\rightarrow 0$, and the leading order term is negative. We remark that the frequency $\omega^{\varepsilon}$ is a real number since the operator in the defect problem can be considered as a real compact perturbation of a self-adjoint periodic operator [@figotin2; @MMMP4].

Therefore, we have proved the following theorem which is the main result of this paper.

Assume that $\delta$ is small enough and the pair $(R,{\varepsilon})$ satisfies one of the two assumptions

-   $R$ small enough and ${\varepsilon}<0$ small enough in magnitude (Dilute regime),

-   $R$ close enough to $1/2$ and ${\varepsilon}>0$ small enough (Non-dilute regime).

Then there exists one frequency value $\omega^{\varepsilon}$ such that the problem has a non-trivial solution. Moreover, in both cases, $\omega^{\varepsilon}$ is slightly above $\omega^*$ and we have $$\begin{aligned}
\omega^{\varepsilon}-\omega_* = 
  \exp\left(-\frac{4\pi^2c_\delta\omega^* R^3}{\delta{\varepsilon}\left(R\|\psi_{\alpha^*}\|^2_{L^2({\partial}D)}  - 2\capacity_{D,\alpha^*} \right)} +O\left(\frac{1}{{\varepsilon}\ln\delta} + 1\right)\right),\end{aligned}$$ when ${\varepsilon}$ and $\delta$ goes to zero.

Since $\omega^{\varepsilon}$ is slightly above $\omega^*=\max_{\alpha}\omega^\alpha$, we have that $\omega^{\varepsilon}$ is located in the bandgap region. This means that the corresponding function $u^{\varepsilon}$, the solution to , should be localized around the defect.

We have shown that it is enough to perturb the size of the defect slightly to create the localized mode in the subwavelength case. This is different from the photonic crystal case where the defect should be large enough to ensure the existence of the localized mode.

There are different physical mechanisms at work in the two different regimes. In the dilute regime, the interactions are weak between the resonators. Decreasing the defect bubble size increases the Minnaert resonance of that single resonator [@first], creating a mode inside the bandgap. In the non-dilute regime, the bubbles are closely separated and have strong interactions. Increasing the defect bubble size reduces the separation, which increases the resonance frequency of that system of bubbles [@doublenegative]. As we will numerically verify in Section \[sec-4\], the shift occurs at $R=1/3$, which corresponds to the case where the bubble radius and the bubble separation are equal.

Proof of Lemma \[lem:estim1\]  {#subsec:proofs_estim}
------------------------------

**Computation of $\Phi^*$ and $\Psi^*$**. Recall that $\ker \mathcal{A}^*(\omega^*,\delta)$ is spanned by $\Psi^*$ and $\ker {\big(\mathcal{A}^*(\omega^*,\delta)\big)^\dagger}$ is spanned by $\Phi^*$. Let us write $$\Psi^{*} = \begin{pmatrix}
\psi_1^*\\ \psi_2^*
\end{pmatrix}, \quad \Phi^{*} = \begin{pmatrix}
\varphi_1^*\\ \varphi_2^*
\end{pmatrix},$$ and assume that $\Psi^*$ and $\Phi^*$ are chosen so that $\int_{{\partial}D}\psi_1^* = 2\pi R$ and $\int_{{\partial}D}{\varphi}_2^* = 2\pi R$. Next we consider the kernel of the limiting operator of $\mathcal{A}^*(\omega,\delta)$ when $\omega\rightarrow0 $, $ \delta\rightarrow0$. Recall that ${\mathcal{S}}_D^\omega$ has a $\log$-singularity as $\omega\rightarrow 0$, given by [(\[eq:Sexpansion\])]{}. Define $$\hat{\mathcal{A}}^*(\omega^*)=\begin{pmatrix}
\hat{\mathcal{S}}^{\omega^*}_D & - \mathcal{S}_D^{\alpha^*,0} 
\\
-\frac{1}{2}I + \mathcal{K}_D^* & 0
\end{pmatrix}.$$ It is known that $\ker(-\frac{1}{2}I+{\mathcal{K}}_D^*)$ is one-dimensional (see, for instance, [@Ammari2007_polarizationBook]). We choose an element $\psi_1^{0,*}\in L^2({\partial}D)$ such that $$\begin{aligned}
{\displaystyle}\bigg(-\frac{1}{2}I+{\mathcal{K}}_D^*\bigg)[\psi_1^{0,*}] = 0,
\quad
{\displaystyle}\int_{{\partial}D} \psi_1^{0,*} =2 \pi R.\end{aligned}$$ Then we define $
{\displaystyle}\psi_2^{0,*} = (\mathcal{S}_D^{\alpha^*,0})^{-1}[\hat{{\mathcal{S}}}_D^{\omega^*}\big[\psi_1^{0,*}]\big]
$ and denote $$\Psi^{0,*} = \begin{pmatrix}
\psi_1^{0,*}\\ \psi_2^{0,*}
\end{pmatrix}.$$ It can be shown that $\Psi^* = \Psi^{0,*} + O(\delta\ln\delta)$ as $\delta\rightarrow 0$. To compute $\Phi^*$, we use the well-known fact that the kernel $\ker(-\frac{1}{2}I+{\mathcal{K}}_D)$ is spanned by the constant function $\chi_{{\partial}D}$. Denote $$\Phi^{0,*} = \begin{pmatrix}
0\\ \chi_{{\partial}D}
\end{pmatrix}.$$ Then it can be shown that $\Phi^* =\Phi^{0,*} + O(\delta) $ as $\delta\rightarrow 0$.

Since $D$ is a disk we can compute $\psi_1^{0,*}$ and $\psi_2^{0,*}$ explicitly. Since $\int_{{\partial}D} \psi_1^{0,*}=2\pi R$, we have $\psi_1^{0,*}= \chi_{{\partial}D}$. Moreover, we have $$\begin{aligned}
\hat{{\mathcal{S}}}^{\omega^*}_{D}[\chi_{{\partial}D}] &= {\mathcal{S}}_D[\chi_{{\partial}D}] + 2\pi R \eta_{\omega^*} \chi_{{\partial}D} \nonumber \\ 
&= \left( R\ln R + 2\pi R \eta_{\omega^*} \right) \chi_{{\partial}D} \nonumber \\ 
&= R \ln\omega^*\chi_{{\partial}D} + O(1). \label{eq:Schi}\end{aligned}$$ Define $$\psi_{\alpha^*} = ({\mathcal{S}}_D^{\alpha^*,0})^{-1}[\chi_{{\partial}D}].$$ Then we get $\psi_2^{0,*} =  R\ln\omega^*\psi_{\alpha^*} + O(1)$ as $\delta\rightarrow 0$. Expand $\psi_2^{0,*}$ in the Fourier basis as $$\psi_2 = \sum_{m\in {\mathbb{Z}}} a_m e^{im\theta}.$$ We will need the first coefficient $a_0$. Using [(\[eq:Schi\])]{} it follows that $$\begin{aligned}
a_0 &= \frac{\ln\omega^*}{2\pi}\int_{\partial D} \left( {\mathcal{S}}_{D}^{\alpha^*}\right) ^{-1} [\chi_{{\partial}D}] {\: \mathrm{d}}\sigma + O(1) \\
&= -\frac{\ln\omega^*}{2\pi}\capacity_{D,\alpha^*}  + O(1).\end{aligned}$$ Finally, using the method from [@thinlayer] section 4.3, we can derive the following improved formula for $\Phi^*$ as $\delta\rightarrow 0$: $$\Phi^* = \begin{pmatrix}
0\\ \chi_{{\partial}D}
\end{pmatrix} - \delta \begin{pmatrix}
\psi_{\alpha^*} \\ 0
\end{pmatrix} + O(\delta^2\ln\delta).$$

**Proof of (i)**. We need the low frequency asymptotics of the operator ${\mathcal{A}}^*(\omega^*,\delta)$. We use the following asymptotic expansions of the Hankel function for small arguments: $$\begin{aligned}
\big(H_0^{(1)}\big)'(z) &= \frac{2i}{\pi }\frac{1}{z} - \frac{i}{\pi} z\ln z +O(z),
\\
\big(H_0^{(1)}\big)''(z) &= -\frac{2i}{\pi }\frac{1}{z^2} - \frac{i}{\pi} \ln z +O(1).\end{aligned}$$ Straightforward computations show that, for small $k$, $$\begin{aligned}
\frac{d}{dk}{\mathcal{S}}_{D}^k[{\varphi}] = \frac{1}{2\pi k} \int_{{\partial}D} {\varphi}(y) d\sigma(y) -\frac{k\ln k}{4\pi}   \int_{{\partial}D} |x-y|^2 {\varphi}(y) d\sigma(y) + O(k),\end{aligned}$$ and $$\begin{aligned}
\frac{d}{dk}({\mathcal{K}}_{D}^k)^*[{\varphi}] &= -\frac{k\ln k}{2\pi}   \int_{{\partial}D} \langle x-y,\nu_x\rangle {\varphi}(y) d\sigma(y) + O(k).\end{aligned}$$ Using low-frequency asymptotics of the quasi-periodic Green’s function $\Gamma^{\alpha,k}$ [@lecturenotes], it follows that $$\begin{aligned}
\frac{d}{dk}{\mathcal{S}}_{D}^{\alpha^*,k}[{\varphi}] = O(k), \quad \mathrm{and} \quad \frac{d}{dk}\left({\mathcal{K}}_{D}^{\alpha^*,k}\right)^*[{\varphi}] = O(1).\end{aligned}$$ Using $\omega = O(\sqrt{\delta})$ we find, on the subspace $V_{00}$, $$\left( \frac{d}{d \omega}{\mathcal{A}}^*(\omega^*,\delta) \Psi^*\right)_{V_{00}} = \begin{pmatrix} O(\frac{1}{\sqrt{\delta}}) \\ \frac{d}{dk}({\mathcal{K}}_{D}^k)^*[\chi_{{\partial}D}] + O(\delta\ln\delta)\end{pmatrix}_{V_{00}}.$$ In total, as $\delta\rightarrow 0$ we have $$\begin{aligned}
\big \langle \Phi^*, \frac{d}{d \omega}{\mathcal{A}}^*(\omega^*,\delta) \Psi^* \big \rangle
&= \langle \chi_{{\partial}D},\frac{d}{d\omega}({\mathcal{K}}_{D}^{\omega})^*[\chi_{{\partial}D}]\rangle + O(\sqrt{\delta})
\\    
&= -\frac{\omega^*\ln \omega^*}{2\pi}  \int_{{\partial}D} \int_{{\partial}D} \langle x-y,\nu_x\rangle  d\sigma(y) d\sigma(x) 
=-2\pi{\omega^*\ln \omega^*} R^3 + O(\sqrt{\delta}).\end{aligned}$$ So, (i) is proved.

**Proof of (ii)**. Using equation \[eq:ADd\], we have $$({\mathcal{A}}_{D}^{\varepsilon})_{V_{mn}\rightarrow V_{m'n'}} = \delta_{mn}\delta_{m'n'} \frac{(-i)\pi R}{2} \begin{pmatrix} J_n(\omega R)H_n^{(1)}(\omega R) & -  J_n(\omega R)H_n^{(1)}(\omega R_d) \frac{J_n(\omega R) }{J_n(\omega R_d)}\\\omega J_n'(\omega R) H_n^{(1)}(\omega R) & - \delta \omega J_n(\omega R) \big(H_n^{(1)}\big)'(\omega R_d) \frac{J_n'(\omega R)}{J_n'(\omega R_d) }\end{pmatrix}.$$ Consequently, the operator $(\mathcal{A}_D^{\varepsilon}- \mathcal{A}_D)$ is given by $$(\mathcal{A}_D^{\varepsilon}- \mathcal{A}_D)_{V_{mn}\rightarrow V_{m'n'}} = \delta_{mn}\delta_{m'n'} \frac{(-i)\pi RJ_n(\omega R)}{2} \begin{pmatrix} 0 & H_n^{(1)}(\omega R) - \frac{J_n(\omega R)H_n^{(1)}(\omega R_d)  }{J_n(\omega R_d)}\\0 & \delta \omega \left(\big(H_n^{(1)}\big)'(\omega R) - \frac{J_n'(\omega R)\big(H_n^{(1)}\big)'(\omega R_d)}{J_n'(\omega R_d)} \right)\end{pmatrix}.$$ Using asymptotic expansions of the Bessel function $J_n(z)$ and the Hankel function $H_n^{(1)}(z)$, for small $z$, straightforward computations show that $$\begin{aligned}
\frac{(-i)\pi R}{2}\frac{J_n(\omega R)}{J_n(\omega R_d)}\left(H_n^{(1)}(\omega R){J_n(\omega R_d)} - J_n(\omega R)H_n^{(1)}(\omega R_d) \right) &= \begin{cases} R\ln\frac{R}{R_d} + O(\omega\ln\omega) \quad &n = 0 \\ -\frac{R}{2n}\left(1-\frac{R^{2n}}{R_d^{2n}}\right) + O(\omega) &n\neq 0\end{cases} \\
&= -{\varepsilon}+ O({\varepsilon}^2 + \omega\ln\omega),\end{aligned}$$ as ${\varepsilon}\rightarrow 0$. Moreover, we have $$\begin{aligned}
\frac{(-i)\pi RJ_n(\omega R)}{2} \delta \omega \left(\big(H_n^{(1)}\big)'(\omega R) - \frac{J_n'(\omega R)\big(H_n^{(1)}\big)'(\omega R_d)}{J_n'(\omega R_d)} \right) &= \begin{cases} \delta \left(1-\frac{R^2}{R_d^2} \right) + O(\delta\omega^2\ln\omega) \quad &n = 0 \\ O(\delta) &n\neq 0\end{cases} \\ 
&= \begin{cases} 2\delta\frac{{\varepsilon}}{R} + O(\delta{\varepsilon}+ \delta\omega^2\ln\omega) \quad &n = 0 \\ O(\delta) &n\neq 0\end{cases}\end{aligned}$$

We are now ready to compute $\left\langle\Phi^*,\big({\mathcal{A}}_{D}^{\varepsilon}(\omega^*,\delta)-{\mathcal{A}}_{D}(\omega^*,\delta)\big){ \Psi^*} \right\rangle$. Using the expressions for $\Phi^*$ and $\Psi^*$, and using $(\omega^*)^2 = O(\delta)$, we have that as $\delta\rightarrow 0$ and ${\varepsilon}\rightarrow 0$, $$\begin{aligned}
\left\langle\Phi^*,\big({\mathcal{A}}_{D}^{\varepsilon}(\omega^*,\delta)-{\mathcal{A}}_{D}(\omega^*,\delta)\big){ \Psi^*} \right\rangle &= \big\langle -\delta\psi_{\alpha^*}, -R{\varepsilon}\ln\omega^* \psi_{\alpha^*} \big\rangle + \big\langle \chi_{{\partial}D},2 a_0 \delta\frac{{\varepsilon}}{R} \chi_{{\partial}D}\big\rangle + O(\delta + {\varepsilon}^2\delta\ln\delta),
\\ 
&= \delta{\varepsilon}\ln\omega^*\left(R\|\psi_{\alpha^*}\|^2_{L^2({\partial}D)}  - 2\capacity_{D,\alpha^*} \right) + O({\varepsilon}\delta + {\varepsilon}^2\delta\ln\delta).\end{aligned}$$ It remains to study the sign of the above expression. Define the factor $S = S(R)$ as $$S(R) = \left(R\|\psi_{\alpha^*}\|^2_{L^2({\partial}D)}  - 2\capacity_{D,\alpha^*} \right).$$ We will show that $S(R) < 0$ in the dilute case while $S(R) > 0$ in the non-dilute case.

*Dilute case.* When the bubbles are small, we can compute $S(R)$ explicitly. Assume $R = \eta$ for some small $\eta$. Then we have [@bandgap] $${\mathcal{S}}_{D}^{\alpha,0} = {\mathcal{S}}_D + O(\eta).$$ Using the fact that $D$ is a circle, it is easily verified that $$\|\psi_{\alpha^*}\|^2_{L^2({\partial}D)} = \frac{2\pi}{R(\ln R)^2} + O(\eta), \quad \capacity_{D,\alpha^*} = -\frac{2\pi}{\ln R} + O(\eta).$$ We then have $$\left(R\|\psi_{\alpha^*}\|^2_{L^2({\partial}D)}  - 2\capacity_{D,\alpha^*} \right) = \frac{2\pi}{\ln R}\left(\frac{1}{\ln R}+2\right),$$ which is negative for small $R$.

*Non-dilute case.* We will show that $S(R)$ is positive for $R$ large enough. We decompose $\psi_{\alpha^*}$ as follows: $$\psi_{\alpha^*} = c\chi_{{\partial}D} + {\varphi}_{\alpha_*}, \quad  \int_{{\partial}D} {\varphi}_{\alpha_*}{\: \mathrm{d}}\sigma = 0.$$ Then it is easily verified that $$c = -\frac{\capacity_{D,\alpha^*}}{2\pi R},$$ so that $$\|\psi_{\alpha^*}\|^2_{L^2({\partial}D)} = \frac{\capacity_{D,\alpha^*}^2}{2\pi R} + \|{\varphi}_{\alpha^*}\|^2_{L^2({\partial}D)} \geq \frac{\capacity_{D,\alpha^*}^2}{2\pi R}.$$ Consequently, $$S(R) \geq \frac{\capacity_{D,\alpha^*}}{2\pi}\left(\capacity_{D,\alpha^*} - 4\pi \right).$$ Moreover, we have that $\capacity_{D,\alpha^*} \rightarrow \infty$ as $R\rightarrow 1/2$. Indeed, using the characterisation found in [@highfrequency], we have that $$\capacity_{D,\alpha^*} = \int_{Y\setminus D} |\nabla v_0|^2 {\: \mathrm{d}}x,$$ where $v_0$ is the harmonic function in $Y\setminus D$ with boundary values $$v_0 = \begin{cases} 1, \quad &\text{on } {\partial}D \\ 0, & \text{on } {\partial}Y.\end{cases}$$ It is easily verified that this integral is unbounded as $R\rightarrow 1/2$. Hence $S(R)$ is positive for $R$ close to $1/2$.

Numerical illustration {#sec-4}
======================

Here we provide numerical examples showing the existence of the subwavelength localized modes.

**Asymptotic formula.** Observe that the $\delta$-error terms in Lemma \[lem:estim1\] only differ from the leading-order terms by a factor of $\ln\delta$, making these expressions unsuitable for numerical computations. Using the same arguments, it is straight-forward to derive the following more precise asymptotic expressions: $$\big\langle \Phi^*, \frac{d}{d \omega}{\mathcal{A}}^*(\omega^*,\delta) \Psi^* \big\rangle = -2\pi{\omega^*R^3 \ln (\omega^*R)} - \pi\omega^*R^3 -\frac{\delta R}{\omega^*}\capacity_{D,\alpha^*} + O(\delta\ln\delta),$$ and $$\left\langle\Phi^*,\big({\mathcal{A}}_{D}^{\varepsilon}(\omega^*,\delta)-{\mathcal{A}}_{D}(\omega^*,\delta)\big){ \Psi^*} \right\rangle = \delta{\varepsilon}\left(\ln R + 2\pi\eta_{\omega^*}\right) \left(R\|\psi_{\alpha^*}\|^2_{L^2({\partial}D)}  - 2\capacity_{D,\alpha^*} \right) + O({\varepsilon}\delta^{3/2}\ln\delta + {\varepsilon}^2\delta\ln\delta).$$

**Implementation.** To verify the asymptotic formula, we also compute the localised frequency by discretizing the operator $\mathcal{M}^{\varepsilon}$ appearing in . Recall that it is given by $$\mathcal{M}^{\varepsilon}(\omega)=\bigg(I+\big({\mathcal{A}}_{D}^{\varepsilon}(\omega,\delta) - {\mathcal{A}}_{D}(\omega,\delta)\big)\frac{1}{(2\pi)^2}\int_{BZ}{\mathcal{A}}^\alpha(\omega,\delta)^{-1} d\alpha\bigg).$$ We use the Fourier basis for the discretisation as in [@bandgap]. Specifically, we use $e^{i n\theta}, n=-N,\ldots,N$ as basis where $N$ is the truncation order. If we increase $N$, then we should get more accurate results.

The operators ${\mathcal{A}}_D$ and ${\mathcal{A}}_D^{\varepsilon}$ are already represented in the Fourier basis in and , respectively. We next consider the operator ${\mathcal{A}}^\alpha$. In [@bandgap], it was shown that the operators $\mathcal{S}^{\alpha,k}$ and ${\partial}\mathcal{S}_D^{\alpha,k}/\partial \nu \big|^\pm_{{\partial}D}$ have the following matrix representations in the Fourier basis: $${\displaystyle}\mathcal{S}_D^{\alpha,k}|_{{\partial}D} \approx (S^\alpha_{m,n})_{m,n=-N}^N, 
\quad 
\frac{{\partial}\mathcal{S}_D^{\alpha,k}}{{\partial}\nu} \big|^\pm_{{\partial}D} \approx
(\widetilde{S}^{\alpha,\pm}_{m,n})_{m,n=-N}^N,$$ where the matrix elements $S_{m,n}^\alpha$ and $\widetilde{S}^{\alpha}_{m,n}$ are given by $$\label{matrixa} \begin{array}{lll}
S^\alpha_{m,n} &=& {\displaystyle}-\frac{i\pi R}{2} J_n (k R)H_n^{(1)}(k R) \delta_{mn} - \frac{i\pi R}{2} J_n(k R)(-1)^{n-m} Q_{n-m}J_m(k R),
\\
{\noalign{\smallskip}}\widetilde{S}^{\alpha,\pm}_{m,n} &=&{\displaystyle}\pm \frac{1}{2}- k \frac{i\pi R}{2}\Big(J_n \cdot (H_n^{(1)})'+J_n'\cdot H_n^{(1)}\Big)(k R) \delta_{mn}- \frac{i\pi R}{2}J_n(k R)(-1)^{n-m} Q_{n-m}k J_m'(k R).
\end{array}$$ Here, $Q_n$ is so called the lattice sum defined by $$Q_n := \sum_{m\in \mathbb{Z}^2, m\neq 0} H_n^{(1)}(k|m|)e^{i n \arg(m)}e^{i m\cdot \alpha}.$$ For an efficient method for computing the lattice sum $Q_n$, we refer the reader to [@Linton2010]. From (\[matrixa\]), the matrix representation of $\mathcal{A}^\alpha$ immediately follows. We used $N=7$ for the truncation order of the Fourier basis.

It remains to consider the integral $\int_{BZ} (\mathcal{A}^\alpha)^{-1} d\alpha$. We compute the two-dimensional integration with respect to $\alpha$ using Gauss quadrature. Then, we finally get the discretization $(M^{\varepsilon}_{m,n})$ of the operator $\mathcal{M}^{\varepsilon}(\omega)$. To find the resonant frequency for the localized mode, we apply Muller’s method [@lecturenotes] to the determinant of the matrix $(M_{m,n}^{\varepsilon})$. The maximum frequency $\omega^*$ of the first band can be used as a good initial guess. The band structure of the unperturbed crystal is computed exactly as in [@bandgap].

For the material parameters, we use $\rho_w=\kappa_w=5000$ and $\rho_b=\kappa_b=1$. In this case, we have $\delta=0.0002$. Moreover, in both examples below we consider the two cases ${\varepsilon}= -0.03R$ and ${\varepsilon}= 0.03R$.

![image](sign.eps){height="5cm"}

![image](dilutefull.eps){height="5.0cm"} ![image](dilutezoom.eps){height="5.0cm"}

In Figure \[fig:SR\], we show the function $S(R)$ over the range $R\in (0,0.4)$. As it can be seen, the function has a zero at $R_0 = 0.326\approx 1/3$. For $R$ below this zero, the crystal is dilute and we expect that ${\varepsilon}$ must be negative in order to have a defect frequency in the bandgap. For $R$ above this zero, we expect that ${\varepsilon}$ must be positive.

**Example 1.** We first consider the dilute case with $R=0.05$. In Figure \[fig:bandgap\_dilute\], we show the computed band structure for this case and the frequency $\omega^{\varepsilon}$ for the localized mode in the case ${\varepsilon}= -0.03R$. In the case ${\varepsilon}=0.03R$, no singular values of $\mathcal{M}^{\varepsilon}(\omega)$ were found. The black circles represents the band structure of the unperturbed crystal ${\mathcal{C}}$. The horizontal red line corresponds to the frequency $\omega^{\varepsilon}$ for the localized mode. The points $\Gamma, X$ and $M$ represent $\alpha=(0,0)$, $\alpha=(\pi,0)$, and $\alpha=(\pi,\pi)$, respectively.

In Figure \[fig:bandgap\_dilute\] (right), we plot the magnified subwavelength band. In the first band, the maximum of $\omega$ is attained at the point $M$ (or $\alpha=\alpha^*=(\pi,\pi)$). The maximum frequency $\omega^*$ of the first band is $\omega^*\approx0.2591$. The defect frequency computed by discretizing the operator is $\omega^{\varepsilon}_d\approx 0.2592$, while the defect frequency computed using the asymptotic formula is $\omega^{\varepsilon}_a\approx 0.2604$. Observe that the defect frequency is exponentially close to $\omega^*$, leading to a high relative error in the computations. However, the frequency $\omega^{\varepsilon}$ is clearly located in the bandgap and in the subwavelength regime. Therefore, the numerical results indicate the existence of a subwavelength localized mode when ${\varepsilon}<0$ but not when ${\varepsilon}>0$.

**Example 2.** We next consider the non-dilute case with $R=0.45$. In Figure \[fig:bandgap\_nondilute\], we show the computed band structure for this non-dilute case, and the frequency $\omega^{\varepsilon}$ for the localized mode in the case ${\varepsilon}= 0.03R$. No defect frequency was found for ${\varepsilon}= -0.03R$ The maximum frequency $\omega^*$ of the first band is $\omega^*\approx0.10806$. The defect frequency computed by discretizing the operator is $\omega^{\varepsilon}_d\approx 0.10807$, while the defect frequency computed using the asymptotic formula is $\omega^{\varepsilon}_a\approx 0.10844$. Again, the frequency $\omega^{\varepsilon}$ is located above the subwavelength band.

Concluding remarks {#sec-5}
==================

In this paper, we have proved for the first time the possibility of localizing waves at the deep subwavelength scale. We have considered a bubbly crystal and produced a localized mode by perturbing the size of one bubble. Moreover, we have proven that the sign of the perturbation depends of the crystal: in the dilute regime the defect bubble should be smaller than the unperturbed bubbles, while in the non-dilute regime the defect bubble should be larger. We have illustrated our findings with numerical experiments. Our results in this paper shed light on the mechanism behind the control and guiding of waves at deep subwavelength scales. In forthcoming works, we plan to investigate the robustness of such subwavelength defect modes with respect to spatial disorder. We will also consider the possibility of guiding waves at deep subwavelength scales using line defects. Another challenging problem we plan to tackle is to prove or disprove the existence of fundamental limits on the Purcell factor in subwavelength bandgap materials.

![image](ndilutefull.eps){height="5.0cm"} ![image](ndilutezoom.eps){height="5.0cm"}

[^1]: Department of Mathematics, ETH Zürich, Rämistrasse 101, CH-8092 Zürich, Switzerland (habib.ammari@math.ethz.ch, brian.fitzpatrick@sam.math.ethz.ch, sanghyeon.yu@sam.math.ethz.ch).

[^2]: Engineering physics, Uppsala University, Sweden (Erik.Orvehed$\_$Hiltunen.0207@student.uu.se).
