Subwavelength phononic bandgap opening in bubbly media

TL;DR

This paper demonstrates the existence of subwavelength phononic bandgaps in bubbly media through analytical and numerical methods, introducing a novel formula for bubble resonance frequencies.

Contribution

It provides an original formula for Minnaert resonance frequencies of arbitrarily shaped bubbles, advancing understanding of phononic bandgap phenomena in bubbly media.

Findings

Existence of subwavelength phononic bandgap in bubbly media

New analytical formula for bubble resonance frequencies

Numerical experiments confirming theoretical predictions

Abstract

The aim of this paper is to show both analytically and numerically the existence of a subwavelength phononic bandgap in bubble phononic crystals. The key is an original formula for the quasi-periodic Minnaert resonance frequencies of an arbitrarily shaped bubble. The main findings in this paper are illustrated with a variety of numerical experiments.

Equations

$$\left\{\begin{array}[]{ll}&\displaystyle\nabla\cdot\frac{1}{\rho}\nabla u+\frac{\omega^{2}}{\kappa}u=0\quad\text{in}\quad Y\backslash D,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr&\displaystyle\nabla\cdot\frac{1}{\rho_{b}}\nabla u+\frac{\omega^{2}}{\kappa_{b}}u=0\quad\text{in}\quad D,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr&\displaystyle u_{+}-u_{-}=0\quad\text{on}\quad\partial D,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr&\displaystyle\frac{1}{\rho}\frac{\partial u}{\partial{\mathbf{\nu}}}\bigg{|}_{+}-\frac{1}{\rho_{b}}\frac{\partial u}{\partial{\mathbf{\nu}}}\bigg{|}_{-}=0\quad\text{on}\quad\partial D,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr&e^{-i\alpha\cdot x}u\,\,\,\mbox{is periodic.}\end{array}\right.$$

$$v=\sqrt{\frac{\kappa }{\rho }},v_b=\sqrt{\frac{{\kappa }_b}{{\rho }_b}},k=\frac{\omega }{v}\text{and}{k}_b=\frac{\omega }{v_b}$$

$$\delta =\frac{{\rho }_b}{\rho }\text{and}\tau =\frac{k_b}{k}=\frac{v}{v_b}=\sqrt{\frac{{\rho }_b\kappa }{\rho {\kappa }_b}}.$$

$$\delta \ll 1,\tau =O(1).$$

$$0\le {\omega }_1^{\alpha }\le {\omega }_2^{\alpha }\le \cdots$$

$$[0,\underset{\alpha }{\max }{\omega }_1^{\alpha }]\cup [\underset{\alpha }{\min }{\omega }_2^{\alpha },\underset{\alpha }{\max }{\omega }_2^{\alpha }]\cup [\underset{\alpha }{\min }{\omega }_3^{\alpha },\underset{\alpha }{\max }{\omega }_3^{\alpha }]\cup \cdots .$$

$${\mathcal{S}}_D^k[\psi ](\mathbf{x}):={\int }_{\partial D}{G}^k(\mathbf{x},\mathbf{y})\psi (\mathbf{y})\mathrm{d}\sigma (\mathbf{y}),$$

$$G^k(\mathbf{x},\mathbf{y}):=-\frac{{\mathrm{e}}^{\mathrm{i}k|\mathbf{x}-\mathbf{y}|}}{4\pi |\mathbf{x}-\mathbf{y}|},$$

$$({\mathcal{K}}_D^k{)}^{\ast }[\psi ](\mathbf{x}):=\text{mbox}p.v.{\int }_{\partial D}\frac{\partial G_k(\mathbf{x},\mathbf{y})}{\partial \nu (\mathbf{x})}\psi (\mathbf{y})\mathrm{d}\sigma (\mathbf{y}),\forall \mathbf{x}\in \partial D.$$

$$\frac{\partial u}{\partial\nu}\Big{|}_{\pm}({\mathbf{x}})=\lim_{t\rightarrow 0^{+}}\langle\nabla u({\mathbf{x}}\pm t\nu({\mathbf{x}})),\nu({\mathbf{x}})\rangle,$$

$$\frac{\partial}{\partial\nu}\Big{|}_{\pm}\mathcal{S}_{D}^{k}[\phi]({\mathbf{x}})=\left(\pm\frac{1}{2}I+(\mathcal{K}_{D}^{k})^{*}\right)[\phi]({\mathbf{x}}),\quad\mbox{a.e.}\leavevmode\nobreak\ {\mathbf{x}}\in\partial D.$$

$$({△}_{\mathbf{x}}+k^2)G^{\alpha ,k}(\mathbf{x},\mathbf{y})=\sum _{n\in {\mathbb{R}}^3}\delta (\mathbf{x}-\mathbf{y}-n)e^{\mathrm{i}n\cdot \alpha },$$

$$G^{\alpha ,k}(\mathbf{x},\mathbf{y})=\sum _{n\in {\mathbb{Z}}^3}\frac{e^{\mathrm{i}(2\pi n+\alpha )\cdot (\mathbf{x}-\mathbf{y})}}{k^2-|2\pi n+\alpha {|}^2},$$

$${\mathcal{S}}_D^{\alpha ,k}[\varphi ](\mathbf{x})={\int }_{\partial D}{G}^{\alpha ,k}(\mathbf{x},\mathbf{y})\varphi (\mathbf{y})d\sigma (\mathbf{y}),\mathbf{x}\in {\mathbb{R}}^3.$$

$$\frac{\partial}{\partial\nu}\Big{|}_{\pm}\mathcal{S}_{D}^{\alpha,k}[\phi]({\mathbf{x}})=\left(\pm\frac{1}{2}I+(\mathcal{K}_{D}^{-\alpha,k})^{*}\right)[\phi]({\mathbf{x}}),\quad\mbox{a.e.}\leavevmode\nobreak\ {\mathbf{x}}\in\partial D,$$

$$({\mathcal{K}}_D^{-\alpha ,k}{)}^{\ast }[\varphi ](\mathbf{x})=\text{mbox}p.v.{\int }_{\partial D}\frac{\partial }{\partial \nu (x)}{G}^{\alpha ,k}(\mathbf{x},\mathbf{y})\varphi (\mathbf{y})d\sigma (\mathbf{y}).$$

$$u=\{\begin{array}{ll}{\mathcal{S}}_D^{\alpha ,k}[\psi ]& \text{in}\text{leavevmode}Y\setminus \bar{D},\\ {\mathcal{S}}_D^{k_b}[{\psi }_b]& \text{in}\text{leavevmode}D,\end{array}$$

$$\mathcal{A}(\omega ,\delta )[\Psi ]=0,$$

$$\mathcal{A}(\omega ,\delta )=\left(\begin{array}{cc}{\mathcal{S}}_D^{k_b}& -{\mathcal{S}}_D^{\alpha ,k}\\ -\frac{1}{2}+{\mathcal{K}}_D^{k_b,\ast }& -\delta (\frac{1}{2}+({\mathcal{K}}_D^{-\alpha ,k}{)}^{\ast })\end{array}\right),\Psi =\left(\begin{array}{c}{\psi }_b\\ \psi \end{array}\right).$$

$$0\le {\omega }_1^{\alpha }\le {\omega }_2^{\alpha }\le \cdots .$$

$$\mathcal{A}(\omega ,0)=\left(\begin{array}{cc}{\mathcal{S}}_D^{k_b}& -{\mathcal{S}}_D^{\alpha ,k}\\ -\frac{1}{2}+{\mathcal{K}}_D^{k_b,\ast }& 0\end{array}\right).$$

$${\mathcal{A}}_0:=\mathcal{A}(0,0)=\left(\begin{array}{cc}{\mathcal{S}}_D& -{\mathcal{S}}_D^{\alpha ,0}\\ -\frac{1}{2}+{\mathcal{K}}_D^{\ast }& 0\end{array}\right),$$

$$\left(-\frac{1}{2}I+{\mathcal{K}}_D^{\ast }\right){\psi }_0=0,{\int }_{\partial D}{\psi }_0=1.$$

$${\mathcal{S}}_D[{\psi }_0]=-{\mathrm{Cap}}_D^{-1}{\text{mathds}1}_{\partial D}.$$

$${\Psi }_0=\left(\begin{array}{c}{\psi }_0\\ {\tilde{\psi }}_0\end{array}\right)\text{and}{\Phi }_0=\left(\begin{array}{c}0\\ {\text{mathds}1}_{\partial D}\end{array}\right),$$

$$\mathcal{A}(\omega ,\delta ):={\mathcal{A}}_0+\mathcal{B}(\omega ,\delta )={\mathcal{A}}_0+\omega {\mathcal{A}}_{1,0}+{\omega }^2{\mathcal{A}}_{2,0}+{\omega }^3{\mathcal{A}}_{3,0}+\delta {\mathcal{A}}_{0,1}+\delta {\omega }^2{\mathcal{A}}_{2,1}+O(|\omega {|}^4+|\delta {\omega }^3|)$$

$${\mathcal{A}}_{1,0}=\left(\begin{array}{cc}{v}_b^{-1}{\mathcal{S}}_{D,1}& 0\\ 0& 0\end{array}\right),{\mathcal{A}}_{2,0}=\left(\begin{array}{cc}{v}_b^{-2}{\mathcal{S}}_{D,2}& -v^{-2}{\mathcal{S}}_{D,1}^{\alpha }\\ v_b^{-2}{\mathcal{K}}_{D,2}^{\ast }& 0\end{array}\right),{\mathcal{A}}_{3,0}=\left(\begin{array}{cc}{v}_b^{-3}{\mathcal{S}}_{D,3}& 0\\ v_b^{-3}{\mathcal{K}}_{D,3}^{\ast }& 0\end{array}\right),$$

$${\mathcal{A}}_{0,1}=\left(\begin{array}{cc}0& 0\\ 0& -(\frac{1}{2}+({\mathcal{K}}_D^{-\alpha ,0}{)}^{\ast })\end{array}\right),{\mathcal{A}}_{2,1}=\left(\begin{array}{cc}0& 0\\ 0& -v^{-2}({\mathcal{K}}_{D,1}^{\alpha }{)}^{\ast }\end{array}\right).$$

$$(\widetilde{{\mathcal{A}}_0}-{\mathcal{P}}_0+\mathcal{B})[{\Psi }_0+{\Psi }_1]=0,$$

$$({\Psi }_1,{\Psi }_0)=0$$

Full text

Mathematics Subject Classification (MSC2000): 35R30, 35C20.

Keywords: Minnaert resonance, bubble, photonic band gap, layer potentials, acoustic waves.

1 Introduction

The past decade has witnessed growing interest in the fabrication of artificially engineered materials to effectively control mechanical waves such as sound waves. Phononic crystals which consist of periodic arrangement of components with controlled spatial size and elastic properties are typical examples. When excited by an acoustic or elastic wave, phononic crystals may exhibit band gaps, or ranges of frequencies in which the wave cannot propagate through their bulk and decaying exponentially. The bandgaps in phononic crystals are due to destructive interference mechanisms from Bragg scattering of the waves, and can be quite wide. For bandgaps to occur, the period of the structure (or the lattice constant) must be of the order of the wavelength and the contrast in the material parameters must be large [4, 3, 5, 13, 15, 20]. This limits the use of phononic crystals in applications targeting low frequencies, because phononic crystals would require impractically large geometries [23, 25].

Based on the realization that composites with locally resonant microstructures can exhibit effective negative elastic parameters at certain frequency ranges, a class of phononic crystals that exhibits bandgaps with lattice constants two orders of magnitude smaller than the wavelength have been fabricated [21]. By varying the size and geometry of the microstructure, it was experimentally shown in [21] (and analytically verified using a simple model in [22]) that one can tune the frequency ranges over which the effective elastic parameters are negative. More recently, resonance has been shown both experimentally and numerically to be another way to prevent elastic waves from propagating in the material in [23].

In this paper, to demonstrate the opening of a subwavelength phononic bandgap, we consider a periodic arrangement of bubbles and exploit their Minnaert resonance [24]. The Minnaert resonant frequency depends on the bulk modulus of the air, the density of the water and the shape of the bubble [24, 17, 1]. In the dilute regime, it has been shown in [10] that around the Minnaert resonant frequency, an effective medium theory can be derived. Furthermore, above the Minnaert resonant frequency, the real part of the effective modulus is negative and consequently, the bubbly fluid behaves as a diffusive media for the acoustic waves. Meanwhile, below the Minnaert resonant frequency, with an appropriate bubble volume fraction, a high contrast effective medium can be obtained, making the superfocusing of waves achievable [9]. These show that the bubbly fluid functions like an acoustic metamaterial and indicate that a subwavelength bandgap opening occurs at the Minneaert resonant frequency [16] . We remark that such behavior is rather analogous to the coupling of electromagnetic waves with plasmonic nanoparticles, which results in effective negative or high contrast dielectric constants for frequencies near the plasmonic resonance frequencies [7, 8].

In this paper, we provide a mathematical and numerical framework for analyzing bandgap opening in bubble phononic crystal at low-frequencies. Through the application of layer potential techniques, Floquet theory, and Gohberg-Sigal theory we derive an original formula for the quasi-periodic Minnaert resonance frequencies of an arbitrarily shaped bubble, along with proving the existence of a subwavelength bandgap and estimating its width. 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 [16]. They pave the mathematical foundation for the analysis of complex-bubble-based phononic crystals that could have more than one structural period and bubbles of different sizes and shapes.

The paper is organized as follows. In Section 2 we formulate the spectral problem for a bubble phononic crystal and introduce some basic results regarding the quasi-periodic Green’s function. In Section 3 we derive an asymptotic formula in terms of the contrast between the densities of the air inside the bubbles and the fluid outside the bubbles. We prove the existence of a subwavelength bandgap and estimate its width. We also consider the dilute regime where the volume fraction of the bubbles is small. In Section 4 we perform numerical simulations in two dimensions to illustrate the main findings of this paper. We make use of the multipole expansion method to compute the subwavelength bandgap. The paper ends with some concluding remarks. In Appendix A, we collect some useful asymptotic formulas for layer potentials in three dimensions. Derivations of the two-dimensional quasi-periodic Minnaert resonances are given in Appendix B. In Appendix C, we briefly describe the basic ideas behind the multipole expansion method.

2 Problem formulation and preliminaries results

We first describe the bubble phononic crystal under consideration. Assume that the bubbles occupy $\cup_{n\in\mathbb{Z}^{d}}(D+n)$ for a bounded and simply connected domain $D$ with $\partial D\in C^{1,s}$ with $0<s<1$. We denote by $\rho_{b}$ and $\kappa_{b}$ the density and the bulk modulus of the air inside the bubbles, respectively, and by $\rho$ and $\kappa$ the corresponding parameters for the background media. To investigate its phononic gap we consider the following $\alpha-$periodic equation in the unit cell $Y=[-1/2,1/2]^{d}$ for $d=2,3$:

[TABLE]

Here, $\partial/\partial{\mathbf{\nu}}$ denotes the outward normal derivative and $|_{\pm}$ denote the limits from outside and inside $D$.

Let

[TABLE]

be respectively the speed of sound outside and inside the bubbles, and the wavenumber outside and inside the bubbles. We also introduce two dimensionless contrast parameters

[TABLE]

By choosing proper physical units, we may assume that the size of the bubble is of order 1. 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 modulus, that is,

[TABLE]

It is known that (2.1) has nontrivial solution for discrete values of $\omega$ such as (see [4])

[TABLE]

and we have the following band structure of propagating frequencies for the given periodic structure:

[TABLE]

In this paper we investigate whether there is a possibility of bandgap opening in this structure.

To do this, we first collect notations and some results regarding the Green function and the quasi-periodic Green’s function for the Helmholtz equation in three dimensions. We refer to [4] and the references therein for the details.

We introduce the single layer potential $\mathcal{S}_{D}^{k}:L^{2}(\partial D)\rightarrow H^{1}(\partial D),H^{1}_{\mathrm{loc}}(\mathbb{R}^{3})$ associated with $D$ and the wavenumber $k$ defined by, $\forall{\mathbf{x}}\in\mathbb{R}^{3},$

[TABLE]

where

[TABLE]

is the Green function of the Helmholtz equation in $\mathbb{R}^{3}$, subject to the Sommerfeld radiation condition. Here, $L^{2}(\partial D)$ is the space of square integrable functions and $H^{1}(\partial D)$ is the standard Sobolev space.

We also define the boundary integral operator $(\mathcal{K}_{D}^{k})^{*}:L^{2}(\partial D)\rightarrow L^{2}(\partial D)$ by

[TABLE]

Here p.v. stands for the Cauchy principal value. We use the notation \frac{\partial}{\partial\nu}\Big{|}_{\pm} indicating

[TABLE]

with $\nu$ being the outward unit normal vector to $\partial D$. Then the following jump formula holds:

[TABLE]

Let $Y$ be the unit cell $[0,1]^{3}$ in $\mathbb{R}^{3}$. For $\alpha\in[-\pi,\pi[^{3}$, the function $G^{\alpha,k}$ is defined to satisfy

[TABLE]

where $\delta$ is the Dirac delta function and $G^{\alpha,k}$ is $\alpha$-quasi-periodic, i.e., $e^{-{\mathrm{i}}\alpha\cdot{\mathbf{x}}}G^{\alpha,k}({\mathbf{x}},{\mathbf{y}})$ is periodic in ${\mathbf{x}}$ with respect to $Y$. It is known that $G^{\alpha,k}$ can be written as

[TABLE]

if $k\neq|2\pi n+\alpha|$ for any $n\in\mathbb{Z}^{3}$.

Let $D$ be a bounded domain in $\mathbb{R}^{3}$ with a connected Lipschitz boundary satisfying $\overline{D}\subset Y$. We define a quasi-periodic single layer potential $\mathcal{S}_{D}^{\alpha,k}$ by

[TABLE]

Then $\mathcal{S}^{\alpha,k}[\phi]$ is an $\alpha$-quasi-periodic function satisfying the Helmholtz equation $(\triangle+k^{2})u=0$. It satisfies a jump formula:

[TABLE]

where $(\mathcal{K}_{D}^{-\alpha,k})^{*}$ is the operator given by

[TABLE]

We remark that it is known that $\mathcal{S}_{D}^{0},\leavevmode\nobreak\ \mathcal{S}_{D}^{\alpha,0}:L^{2}(\partial D)\rightarrow H^{1}(\partial D)$ are invertible for $\alpha\neq 0$; see [4].

3 Subwavelength bandgaps

We use layer potentials to represent the solution to the scattering problem (2.1). We look for a solution $u$ of (2.1) of the form

[TABLE]

for some surface potentials $\psi,\psi_{b}\in L^{2}(\partial D)$. Using the jump relations for the single layer potentials, one can show that (2.1) is equivalent to the boundary integral equation

[TABLE]

where

[TABLE]

Throughout the paper, we denote by $\mathcal{H}=L^{2}(\partial D)\times L^{2}(\partial D)$ and by $\mathcal{H}_{1}=H^{1}(\partial D)\times L^{2}(\partial D)$, and use $(\cdot,\cdot)$ for the inner product in $L^{2}$ spaces and $\|\cdot\|_{\mathcal{H}}$ for the norm in $\mathcal{H}$. It is clear that $\mathcal{A}(\omega,\delta)$ is a bounded linear operator from $\mathcal{H}$ to $\mathcal{H}_{1}$, i.e. $\mathcal{A}(\omega,\delta)\in\mathcal{B}(\mathcal{H},\mathcal{H}_{1})$. Moreover, we can check that the characteristic values of $\mathcal{A}(\omega,\delta)$ can be written as

[TABLE]

We first look at the limiting case when $\delta=0$. The operator $\mathcal{A}(\omega,\delta)$ is a perturbation of

[TABLE]

We see that $\omega_{0}$ is a characteristic value of $\mathcal{A}(\omega,0)$ if only if $(\omega_{0}/v_{b})^{2}$ is a Neumann eigenvalue of $D$ or $(\omega_{0}/v)^{2}$ is a Dirichlet eigenvalue of $Y\backslash D$ with $\alpha$-quasiperiodicity on $\partial Y$. Since zero is a Neumann eigenvalue of $D$, $\omega=0$ is a characteristic value for the operator-valued analytic function $\mathcal{A}(\omega,0)$. Besides, note that there is a positive lower bound for other Neumann eigenvalues of $D$ and all the Dirichlet eigenvalues of $Y\backslash D$ with $\alpha$-quasiperiodicity on $\partial Y$, we can conclude the following result by the Gohberg-Sigal theory [4, 14].

Lemma 3.1**.**

For any $\delta$ sufficiently small, there exists one and only one characteristic value $\omega_{0}=\omega_{0}(\delta)$ in a neighborhood of the origin in the complex plane to the operator-valued analytic function $\mathcal{A}(\omega,\delta)$. Moreover, $\omega_{0}(0)=0$ and $\omega_{0}$ depends on $\delta$ continuously.

3.1 The asymptotic behavior of $\omega_{1}^{\alpha}$

In this section we assume $\alpha\neq 0$. We define

[TABLE]

Here we set for brevity ${\mathcal{S}}_{D}:={\mathcal{S}}_{D}^{k=0}$, ${\mathcal{K}}_{D}^{*}:={\mathcal{K}}_{D}^{k=0,*}$. We denote by ${\mathds{1}}_{\partial D}\in H^{1}(\partial D)$ the constant function on $\partial D$ with value $1$, and by $\mathcal{A}_{0}^{*}:{\mathcal{H}}_{1}\rightarrow{\mathcal{H}}$ the adjoint of $\mathcal{A}_{0}$. We choose an element $\psi_{0}\in L^{2}(\partial D)$ such that

[TABLE]

We define the capacity of the set $D$, ${\mathrm{Cap}}_{D}$, by

[TABLE]

Then we can easily check that ${\mathrm{Ker}}(\mathcal{A}_{0})$ and ${\mathrm{Ker}}(\mathcal{A}_{0}^{*})$ are spanned respectively by

[TABLE]

where $\tilde{\psi}_{0}=(\mathcal{S}_{D}^{\alpha,0})^{-1}\mathcal{S}_{D}[\psi_{0}]$.

We now perturb $\mathcal{A}_{0}$ by a rank-1 operator $\mathcal{P}_{0}$ from $\mathcal{H}$ to $\mathcal{H}_{1}$ given by $\mathcal{P}_{0}[\Psi]:=(\Psi,\Psi_{0})\Phi_{0},$ and denote it by $\tilde{\mathcal{A}_{0}}=\mathcal{A}_{0}+\mathcal{P}_{0}$.

Lemma 3.2**.**

The followings hold:

1. (i)

$\tilde{\mathcal{A}_{0}}[\Psi_{0}]=\|\Psi_{0}\|^{2}\Phi_{0}$, $\tilde{\mathcal{A}_{0}}^{*}[\Phi_{0}]=\|\Phi_{0}\|^{2}\Psi_{0}$. 2. (ii)

The operator $\tilde{\mathcal{A}_{0}}$ and its adjoint $\tilde{\mathcal{A}_{0}}^{*}$ are invertible in $\mathcal{B}(\mathcal{H},\mathcal{H}_{1})$ and $\mathcal{B}(\mathcal{H}_{1},\mathcal{H})$, respectively.

Proof.

By construction, and the fact that $\mathcal{S}_{D}$ is bijective from $L^{2}(\partial D)$ to $H^{1}(\partial D)$ [2], we can show that $\tilde{\mathcal{A}_{0}}$ (hence $\tilde{\mathcal{A}_{0}}^{*}$) is bijective. The fact that $\tilde{{\mathcal{A}}_{0}}[\Psi_{0}]=\|\Psi_{0}\|^{2}\Phi_{0}$ is direct. Finally, by noticing that $\mathcal{P}_{0}^{*}[\theta]=(\theta,\Phi_{0})\Psi_{0}$, it follows that $\tilde{\mathcal{A}_{0}}^{*}[\Phi_{0}]=\mathcal{P}_{0}^{*}[\Phi_{0}]=\|\Phi_{0}\|^{2}\Psi_{0}.$ ∎

Using the results in Appendix A, we can expand $\mathcal{A}(\omega,\delta)$ as

[TABLE]

where

[TABLE]

[TABLE]

Since $\tilde{\mathcal{A}_{0}}=\mathcal{A}_{0}+\mathcal{P}_{0}$, the equation (3.2) is equivalent to

[TABLE]

where

[TABLE]

Observe that the operator $\tilde{\mathcal{A}_{0}}+\mathcal{B}$ is invertible for sufficiently small $\delta$ and $\omega$. Applying $(\tilde{\mathcal{A}_{0}}+\mathcal{B})^{-1}$ to both sides of the above equation leads to

[TABLE]

Using the condition $(\Psi_{1},\Psi_{0})=0$, we deduce that (3.2) has a nontrivial solution if and only if

[TABLE]

Let us calculate $A(\omega,\delta):=\widetilde{A}(\omega,\delta)\left\|\Phi_{0}\right\|^{2}$. Using the Neumann series

[TABLE]

and the fact that $\tilde{\mathcal{A}_{0}}^{-1}[\Phi_{0}]=\left\|\Psi_{0}\right\|^{-2}\Psi_{0}$ and $(\tilde{\mathcal{A}_{0}}^{*})^{-1}[\Psi_{0}]=\left\|\Phi_{0}\right\|^{-2}\Phi_{0}$, we obtain that

[TABLE]

It is clear that $\mathcal{A}_{1,0}^{*}[\Phi_{0}]=0$. Consequently, the expression simplifies into

[TABLE]

We now calculate the five remaining terms.

$\bullet$** Calculation of $\left(\mathcal{A}_{2,0}[\Psi_{0}],\Phi_{0}\right)$**. Using the first point of Lemma A.1, we get

[TABLE]

where we used the fact that ${\mathcal{S}}_{D}[\psi_{0}]({\mathbf{x}})=-{\mathrm{Cap}}_{D}^{-1}$ for all ${\mathbf{x}}\in D$.

$\bullet$** Calculation of $\left(\mathcal{A}_{3,0}[\Psi_{0}],\Phi_{0}\right)$**. Similarly, using the second point of Lemma A.1, we get

[TABLE]

$\bullet$** Calculation of $\left(\mathcal{A}_{0,1}[\Psi_{0}],\Phi_{0}\right)$**. We directly have

[TABLE]

$\bullet$** Calculation of $\left(\mathcal{A}_{0,1}\tilde{\mathcal{A}_{0}}^{-1}\mathcal{A}_{1,0}[\Psi_{0}],\Phi_{0}\right)$**. We have

[TABLE]

Let us calculate $\tilde{\mathcal{A}}_{0}^{-1}\begin{pmatrix}{\mathds{1}}_{\partial D}\\ 0\end{pmatrix}$. We look for $(a\psi_{0},b\tilde{\psi}_{0})\in{\mathcal{H}}$ so that

[TABLE]

By solving the above equations directly, we obtain

[TABLE]

It follows that

[TABLE]

$\bullet$** Calculation of $\left(\mathcal{A}_{2,0}\tilde{\mathcal{A}_{0}}^{-1}\mathcal{A}_{1,0}[\Psi_{0}],\Phi_{0}\right)$**. Using similar calculations, we obtain

[TABLE]

$\bullet$** Conclusion.** Considering the above the results, we can derive from (3.9) that

[TABLE]

where

[TABLE]

and

[TABLE]

We now solve $A(\omega,\delta)=0$. It is clear that $\delta=O(\omega^{2})$ and thus $\omega_{0}(\delta)=O(\sqrt{\delta})$. We write $\omega_{0}(\delta)=a_{1}\delta^{\frac{1}{2}}+a_{2}\delta+O(\delta^{\frac{3}{2}})$, and get

[TABLE]

From the coefficients of the $\delta$ and $\delta^{\frac{3}{2}}$ terms, we obtain

[TABLE]

which yields

[TABLE]

Therefore, we obtain

Theorem 3.1**.**

For $\alpha\neq 0$ and sufficiently small $\delta$, we have

[TABLE]

where $\omega_{M}=\sqrt{\frac{\delta v_{b}^{2}{\mathrm{Cap}}_{D}}{|D|}}$ is the (free space) Minnaert resonant frequency.

Let us define the $\alpha$-quasi-periodic capacity by

[TABLE]

Then we have

[TABLE]

and (4) is written as

[TABLE]

with $\omega_{M,\alpha}=\sqrt{\frac{\delta v_{b}^{2}{\mathrm{Cap}}_{D,\alpha}}{|D|}}.$ We can see that

[TABLE]

as $\alpha\rightarrow 0$ because $\left(1/2+(\mathcal{K}_{D}^{-\alpha,0})^{*}\right)(\mathcal{S}_{D}^{\alpha,0})^{-1}[{\mathds{1}}_{\partial D}]\rightarrow 0$ and so ${\mathrm{Cap}}_{D,\alpha}\rightarrow 0$ as $\alpha\rightarrow 0$.

We define $\omega_{1}^{*}:=\max_{\alpha}\omega_{M,\alpha}$. Then we deduce the following regarding a bandgap opening.

Theorem 3.2**.**

For every $\epsilon>0$, there exists $\delta_{0}>0$ and $\tilde{\omega}>\omega_{1}^{*}+\epsilon$ such that

[TABLE]

for $\delta<\delta_{0}$.

Proof.

Using $\omega_{1}^{0}=0$ and the continuity of $\omega_{1}^{\alpha}$ in $\alpha$ and $\delta$, we get $\alpha_{0}$ and $\delta_{1}$ such that $\omega_{1}^{\alpha}<\omega_{1}^{*}$ for every $|\alpha|<\alpha_{0}$ and $\delta<\delta_{1}$. Following the derivation of (4), we can check that it is valid uniformly in $\alpha$ as far as $|\alpha|\geq\alpha_{0}$. Thus there exists $\delta_{0}<\delta_{1}$ such $\omega_{1}^{\alpha}\leq\omega_{1}^{*}+\epsilon$ for $|\alpha|\geq\alpha_{0}$. We have shown that $\max_{\alpha}\omega_{1}^{\alpha}\leq\omega_{1}^{*}+\epsilon$ for sufficiently small $\delta$. To have $\min_{\alpha}\omega_{2}^{\alpha}>\omega_{1}^{*}+\epsilon$ for small $\delta$, it is enough to check that $\mathcal{A}(\omega,\delta)$ has no small characteristic value other than $\omega_{1}^{\alpha}$. For $\alpha$ away from [math], we can see that it is true following the proof of Theorem 3.1. If $\alpha=0$, we have

[TABLE]

near $\omega_{2}^{0}$ with $\delta=0$. Since $\omega_{2}^{0}\neq 0$, we have $\omega_{2}^{0}(\delta)>\omega_{1}^{*}+\epsilon$ for sufficiently small $\delta$. Finally, using the continuity of $\omega_{2}^{\alpha}$ in $\alpha$, we obtain $\min_{\alpha}\omega_{2}^{\alpha}>\omega_{1}^{*}+\epsilon$ for small $\delta$. This completes the proof. ∎

3.2 Dilute case

We emphasize that our calculations in the previous part hold even for the dilute case as long as $\delta/\eta^{2}$ is small where $\eta$ is the diameter of $D$.

We state an asymptotic behavior of ${\mathrm{Cap}}_{D,\alpha}$ when $D=\eta B$ for a small $\eta$. Note that ${\mathrm{Cap}}_{D}={\mathrm{Cap}}_{\eta B}=\eta{\mathrm{Cap}}_{B}$. Fix $c>0$, the following holds.

Lemma 3.3**.**

For $|\alpha|>c>0$, we have

[TABLE]

where $R_{\alpha}({\mathbf{x}}):=G^{\alpha,0}({\mathbf{x}})-G^{0}({\mathbf{x}})$.

Proof.

Since $R_{\alpha}({\mathbf{x}})$ is smooth and $R_{\alpha}({\mathbf{x}})=R_{\alpha}(0)+O(|{\mathbf{x}}|)$ as $|{\mathbf{x}}|\rightarrow 0$, we have

[TABLE]

with $\tilde{\phi}({\mathbf{x}}):=\phi(\eta{\mathbf{x}})$. Then

[TABLE]

and so

[TABLE]

∎

By this approximation we have $\omega_{M,\alpha}\approx\omega_{M}$ for $\alpha$ away from [math] and so $\omega_{1}^{*}\approx\omega_{M}$. Combined with Theorem 3.2 this means that there is a band gap opening slightly above the Minneart resonance frequency for a single bubble. It is coherent with results in [10] showing an effective medium theory for the bubbly fluid as the number of bubbles tends to infinity. It is shown there that near and above the Minnaert resonant frequency, the obtained effective media can have a negative bulk modulus.

4 Numerical illustrations

Recall the formula for the $\alpha$-quasi-periodic Minneart resonance:

[TABLE]

We want to compare $\omega_{approx}^{\alpha}:=\omega_{M}\sqrt{\frac{{\mathrm{Cap}}_{D,\alpha}}{{\mathrm{Cap}}_{D}}}$ with the true $\alpha$-quasi-periodic resonance $\omega_{exact}^{\alpha}$, which can be obtained through direct calculation of the minimum characteristic value of the operator $\mathcal{A}(\omega,\delta)$ in (3.2) using Muller’s method [1].

We set the density and the bulk modulus of the bubbles to be $\rho_{b}=1$ and $\kappa_{b}=1$, respectively. In order to confirm that the formula becomes accurate in the appropriate regime, which features similar wavenumbers inside and outside the bubbles along with, in particular, a high contrast in the bulk modulii, we take the density and the bulk modulus of the background material to be $\rho=\kappa=\delta^{-1}\in(10,1000)$. We assume that the bubble represented by $D$ is a disk of radius $R=0.0125$. In Figure 3 we plot $\omega_{approx}^{\alpha}$ and $\omega_{exact}^{\alpha}$ against the contrast $\delta^{-1}$ and it is clear that the formula provides a highly accurate approximation when the contrast is sufficiently large.

Next we present numerical examples to illustrate subwavelength bandgap openings. As $D$ is a disk of radius $R$, we apply the multipole expansion method for computing the band structure (for the details, we refer to Appendix C). As described in [11], the quasi-periodic Green’s function is unsuitable for bandgap calculations due to empty resonance phenomenon. Therefore, we make use of the multiple expansion method which is efficient in the case of disk-shaped bubbles.

We first consider the dilute case. We set $R=0.05$, $\rho=\kappa=5000$ and $\rho_{b}=\kappa_{b}=1$. In this case, we have $\delta=0.0002$. Figure 1 shows the computed band structure. The points $\Gamma,X$ and $M$ represent $\alpha=(0,0)$, $\alpha=(\pi,0)$ and $\alpha=(\pi,\pi)$, respectively. We plot the first two characteristic values $\mathcal{A}(\omega,\delta)$ along the boundary of the triangle $\Gamma XM$. It can be seen that a subwavelength bandgap in the spectrum of $\mathcal{A}(\omega,\delta)$ does exist. Moreover, the bandgap between the first two bands is quite large. It is also worth mentioning that, by zooming the subwavelength bandgap (on the right in Figure 1), one can see that $\omega_{1}^{\ast}$ is attained at the point $M$ (that is, $\alpha=(\pi,\pi)$). We used $N=7$ for the truncation order of cylindrical waves. Further numerical experiments indicate that this phenomenon is independent of the bubble radius or position.

Next we consider a non-dilute regime. We set $R=0.25$ and $\rho=\kappa=1000$ and $\rho_{b}=\kappa_{b}=1$. In this case we have $\delta=0.001$. Figure 2 shows the computed band structure. Again, a subwavelength bandgap can be observed. We used $N=3$ for the truncation order of cylindrical waves in the multipole expansion method.

Finally, in order to verify our conclusion from Lemma 3.3, namely that $\omega_{1}^{\ast}=\max_{\alpha}\omega_{M,\alpha}\approx\omega_{M}$ when $\alpha$ is non-zero, we fix the contrast to be $\delta^{-1}=1000$ and observe $\omega_{1}^{\ast}$ and $\omega_{M}$ over a range of bubble sizes in Figure 4.

5 Concluding remarks

In this paper we have proved the existence of a subwavelength bandgap opening in bubble phononic crystals. We have illustrated our main findings with a variety of numerical experiments. We have also covered recently proved results on the effective medium theory in the dilute case. In a forthcoming work we will use the Bloch wave decomposition for homogenizing [12, 19] the bubble crystal near the maximum over $\alpha$ of $\omega_{1}^{\alpha}$ in the general case. Our aim is to prove that above such a frequency the crystal behaves like a material with a negative density while below it it behaves like a high contrast material, leading to superfocusing of acoustic waves.

Appendix A Some asymptotic expansions

We recall some basic asymptotic expansions for the layer potentials in three dimensions from [4] (see also the appendix in [6]).

A.1 Asymptotic expansions of $\mathcal{S}_{D}^{k}$ and $\mathcal{K}_{D}^{k,*}$

We expand the Green’s function $G_{k}$ with

[TABLE]

In particular, $G_{1}({\mathbf{x}})=-\frac{{\mathrm{i}}}{4\pi}$. Developing in power of $k$ the equation $(\Delta+k^{2})G_{k}=\delta_{0}$ leads to

[TABLE]

From (A.1), we decompose the single layer potential as

[TABLE]

where the convergence holds in $\mathcal{B}(L^{2}(\partial D),H^{1}(\partial D))$. Similarly, the asymptotic expansion for the operator $\mathcal{K}_{D}^{k,*}$ is

[TABLE]

where the convergence holds in $\mathcal{B}(L^{2}(\partial D),L^{2}(\partial D))$. Using (A.2), we deduce the following useful identities.

Lemma A.1**.**

It holds:

[TABLE]

A.2 Asymptotic expansions of $\mathcal{S}_{D}^{\alpha,k},\leavevmode\nobreak\ (\mathcal{K}_{D}^{\alpha,k})^{*}$

For the $\alpha$-quasi-periodic Green’s function $G^{\alpha,k}$, we have

[TABLE]

when $\alpha\neq 0$, and $k\rightarrow 0$.

From (A.5), we decompose the single layer potential as

[TABLE]

where the convergence holds in $\mathcal{B}(L^{2}(\partial D),H^{1}(\partial D))$. Similarly, the asymptotic expansion for the operator $\left(\mathcal{K}_{D}^{-\alpha,k}\right)^{*}$ is

[TABLE]

where the convergence holds in $\mathcal{B}(L^{2}(\partial D),L^{2}(\partial D))$.

Appendix B The two-dimensional case

The aim of this appendix is to check that formula (4) holds in the two-dimensional case, where $\omega_{M}$ is the (free space) Minnaert resonant frequency and $c_{2}$ is defined by (3.13). Note that for $\alpha\neq 0$, the quasi-periodic single layer operator $\mathcal{S}_{D}^{\alpha,0}:L^{2}(\partial D)\rightarrow H^{1}(\partial D)$ is invertible. Moreover, the definitions (3.5) and (3.15) of both the capacity and the $\alpha$-quasi-periodic capacity remain valid.

Using the asymptotic expansions in [1, Appendix A] as $k\rightarrow 0$,

[TABLE]

where for $\psi\in L^{2}(\partial D)$

[TABLE]

and

[TABLE]

with

[TABLE]

and $\gamma$ being the Euler constant.

Therefore, in the two-dimensional case the asymptotic expansion (3.6) should be replaced with

[TABLE]

where $\mathcal{A}_{0}$ is defined by (3.4),

[TABLE]

and

[TABLE]

Using the definition of the free space Minnaert resonance in dimension two in [1, Theorem B1], it is not difficult to see that (4) holds.

Appendix C Multipole expansion method

When $D$ is a circular disk of radius $R$, the integral equation admits an explicit representation. In this case, the solution can be represented as a sum of cylindrical waves $J_{n}(kr)e^{{\mathrm{i}}n\theta}$ or $H_{n}^{(1)}(kr)e^{{\mathrm{i}}n\theta}$. Here we give a multipole expansion interpretation of the integral operator $\mathcal{A}$, which leads to an efficient numerical scheme for computing its bandgap structure.

Recall that, for each fixed $k,\alpha$, we have to find a characteristic value of $\mathcal{A}(\omega,\delta)$ defined by

[TABLE]

From the above expression, we see that $\mathcal{A}(\omega,\delta)$ is represented in terms of the single layer potential only. So it is enough to derive a multipole expansion version of the single layer potential.

Let us first consider the single layer potential $\mathcal{S}_{D}^{k}[\varphi]$ for a single disk $D$. We adopt the polar coordinates $(r,\theta)$. Then, since $D$ is a circular disk, the density function $\varphi=\varphi(\theta)$ is a $2\pi$-periodic function. It admits the following Fourier series expansion:

[TABLE]

for some coefficients $a_{n}$. Hence we only need to compute $u:=\mathcal{S}_{D}^{k}[e^{{\mathrm{i}}n\theta}]$ which satisfies

[TABLE]

The above equation can be easily solved by the separation of variables technique in polar coordinates. It gives

[TABLE]

where $c=\frac{-{\mathrm{i}}\pi R}{2}$.

Now we compute the quasi-periodic single layer potential $\mathcal{S}_{D}^{\alpha,k}[e^{{\mathrm{i}}n\theta}]$. Since

[TABLE]

we have

[TABLE]

Here, $D+m$ means a translation of the disk $D$ by $m$ and $(r_{m},\theta_{m})$ is the polar coordinates with respect to the center of $D+m$. By applying the following addition theorem:

[TABLE]

we obtain

[TABLE]

where $Q_{n}$ is so called the lattice sum defined by

[TABLE]

So, from (C.3) and (C.4), we finally obtain an explicit representation of $\mathcal{S}_{D}^{\alpha,k}$. For an efficient method for computing the lattice sum $Q_{n}$, see [18].

In the numerical computations, we should consider the truncated series

[TABLE]

instead of $\mathcal{S}_{D}^{\alpha,k}[\varphi]=\sum_{n\in\mathbb{Z}}a_{n}\mathcal{S}_{D}^{\alpha,k}[e^{{\mathrm{i}}n\theta}]$ for some large enough $N\in\mathbb{N}$. Then, using $e^{{\mathrm{i}}n\theta}$ as basis, we have the following matrix representation of the operator $\mathcal{S}^{\alpha,k}$:

[TABLE]

where $S_{m,n}$ is given by

[TABLE]

Similarly, we also have the following matrix representation for $\frac{\partial\mathcal{S}_{D}^{\alpha,k}}{\partial\nu}|_{\partial D}^{\pm}$:

[TABLE]

where $S^{\prime\pm}_{m,n}$ is given by

[TABLE]

The matrix representation of $\mathcal{A}(\omega,\delta)$ immediately follows.