You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
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 (<PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>, <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>, <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>).

[^2]: Engineering physics, Uppsala University, Sweden (Erik.Orvehed$\<PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>).