A High-Frequency Homogenization Approach Near the Dirac Points in Bubbly Honeycomb Crystals

Abstract

In (Ammari et al. in SIAM J Math Anal. arXiv:1811.03905), the existence of a Dirac dispersion cone in a bubbly honeycomb phononic crystal comprised of bubbles of arbitrary shape is shown. The aim of this paper is to prove that, near the Dirac points, the Bloch eigenfunctions is the sum of two eigenmodes. Each eigenmode can be decomposed into two components: one which is slowly varying and satisfies a homogenized equation, while the other is periodic across each elementary crystal cell and is highly oscillating. The slowly oscillating components of the eigenmodes satisfy a system of Dirac equations. Our results in this paper prove for the first time a near-zero effective refractive index near the Dirac points for the plane-wave envelopes of the Bloch eigenfunctions in a sub-wavelength metamaterial. They are illustrated by a variety of numerical examples. We also compare and contrast the behaviour of the Bloch eigenfunctions in the honeycomb crystal with that of their counterparts in a bubbly square crystal, near the corner of the Brillouin zone, where the maximum of the first Bloch eigenvalue is attained.

1 Introduction

Metamaterials are a novel group of materials designed to have special wave characteristics such as bandgaps, negative refractive indices, or sub-wavelength scale resolution in imaging. There have also been demonstrations of materials with near-zero refractive indices. These materials have a wide number of applications, including low-loss bending transmission, invisibility cloaking, and zero phase-shift propagation [12, 14, 21, 29, 37].

The first near-zero refractive index phononic crystal was theoretically demonstrated in [29], where the effective mass density and reciprocal bulk modulus were shown to vanish simultaneously. This near-zero effective refractive index property is a consequence of the existence of a Dirac dispersion cone. The double-zero property is possible because the Dirac cone is located at the centre $\Gamma$ of the Brillouin zone. Single-zero properties have been studied for other locations of the Dirac cone; however, these materials exhibit a low transmittance making them less desirable for applications [14, 19, 20]. At high frequencies, the Bloch eigenmodes will oscillate on the microscale of the metamaterial, suggesting that an “effective parameter” description of the material is overly simplified. Nevertheless, as will be shown in this paper, an effective equation for the envelope of each of these Bloch eigenmodes can be derived.

Metamaterials with Dirac singularities have been experimentally and numerically studied in [33, 35, 36]. Proofs of the existence of a Dirac cone at the symmetry point K in honeycomb lattice structures and mathematical analyses of their properties are provided in [5, 11, 13, 15,16,17, 24, 25, 32, 34]. In [14], time-dependent material parameters are used to move the Dirac cone from the point K to the centre $\Gamma$ of the Brillouin zone, enabling a double-zero refractive index.

Sub-wavelength resonators are the building blocks of metamaterials. In acoustics, a gas bubble in a liquid is known to have a resonance frequency corresponding to wavelengths which are several orders of magnitude larger than the bubble [2, 31]. This opens up the possibility of creating small-scaled acoustic metamaterials known as sub-wavelength metamaterials, whereby the operating frequency corresponds to wavelengths much larger than the device size. The simplicity of the gas bubble makes bubbly media an ideal model for sub-wavelength metamaterials. Many experimentally observed phenomena in bubbly media [22, 23, 26,27,28, 30, 38] have been rigorously explained in [1, 3, 5,6,7, 10]. In particular, in [5], a bubbly honeycomb crystal is considered, and a Dirac dispersion cone centred at the symmetry point K in the Brillouin zone is shown to exist.

In this paper, we prove the near-zero effective index property around the point K in a bubbly honeycomb crystal at the deep sub-wavelength scale. We will develop a homogenization theory that captures both the macroscopic behaviour of the eigenmodes and the oscillations in the microscopic scale, and demonstrate that the near-zero property holds in the macroscopic scale.

In the homogenization theory of metamaterials, the goal is to map the metamaterial to a homogeneous material with some effective parameters. It has previously been demonstrated that this approach does not apply in the case of bubbly crystals at “high” frequencies, i.e. away from the centre $\Gamma$ of the Brillouin zone. In [9], it is shown that around the symmetry point M in the Brillouin zone of a bubbly crystal with a square lattice, the Bloch eigenmodes display oscillatory behaviour on two distinct scales: small scale oscillations on the order of the size of individual bubbles, while simultaneously the plane-wave envelope oscillates at a much larger scale and satisfies a homogenized equation. Analogously, we expect the standard homogenization approach to fail for the honeycomb crystal and seek instead a homogenized equation for the envelopes of the eigenmodes. We will demonstrate that this is a near-zero refractive index homogenized equation near the Dirac points. Moreover, we will compare our results with the case of a square lattice crystal, which does not have a linear dispersion relation around the symmetry points of the Brillouin zone, and consequently cannot have effective near-zero refractive index.

This paper is organized as follows: in Section 2, we present the eigenvalue problem of the bubbly honeycomb crystal, and state the main results of [5]. In Section 3, we use layer-potential techniques to compute the Bloch eigenfunctions close to the point K in the asymptotic limit of high density contrast. In Section 4, we decompose the Bloch eigenfunctions as the sum of two eigenmodes, each with a slowly oscillating plane-wave envelope. We also derive a system of Dirac equations satisfied by the slowly oscillating components of the eigenmodes in the vicinity of the Dirac points, which generalizes to wave propagation in sub-wavelength resonant structures the result obtained in [18] for the Schrödinger equation. The main result of this paper is stated in Theorem 4.3, where the Bloch eigenmodes are shown to exhibit the two-scale behaviour as described above. In Section 5, we numerically illustrate Theorem 4.3. We show that the macroscopic plane-wave envelope in the honeycomb crystal has a lower order of oscillations compared to the Bloch eigenfunction of the square crystal. Moreover, we demonstrate the near-zero effective refractive index property of the honeycomb crystal, and compare the behaviour of the Bloch eigenfunctions with those of a square crystal. Finally, in Section 6, we summarise the main results of this paper and briefly discuss the remaining challenges in the field.

2 Problem Statement and Preliminaries

In this section, we describe the honeycomb lattice and state the main results of [5].

2.1 Problem Formulation

We consider a two-dimensional infinite honeycomb crystal in two dimensions depicted in Fig. 1. Define the lattice $\Lambda$ generated by the lattice vectors

$$\begin{aligned} l_1 = L\left( \frac{\sqrt{3}}{2}, \frac{1}{2} \right) ,~~l_2 = L\left( \frac{\sqrt{3}}{2}, -\frac{1}{2}\right) , \end{aligned}$$

where L is the lattice constant. Denote by Y a fundamental domain of the given lattice. Here, we take

$$\begin{aligned} Y:= \left\{ s l_1+ t l_2 ~|~ 0 \le s,t \le 1 \right\} . \end{aligned}$$

Define the three points $x_0, x_1,$ and $x_2$ as

$$\begin{aligned} x_0 = \frac{l_1 + l_2}{2}, \quad x_1 = \frac{l_1+l_2}{3}, \quad x_2 = \frac{2(l_1 + l_2)}{3}. \end{aligned}$$

Fig. 1

Illustration of the bubbly honeycomb crystal and quantities in the fundamental domain Y

We will consider a general shape of the bubbles, under certain symmetry assumptions. Let $R_0$ be the rotation around $x_0$ by $\pi$, and let $R_1$ and $R_2$ be the rotations by $-\frac{2\pi }{3}$ around $x_1$ and $x_2$, respectively. These rotations can be written as

$$\begin{aligned} R_1 x = Rx+l_1, \quad R_2 x = Rx + 2l_1, \quad R_0 x = 2x_0 - x, \end{aligned}$$

where R is the rotation by $-\frac{2\pi }{3}$ around the origin. Moreover, let $R_3$ be the reflection across the line $p = x_0 + \mathbb {R}e_2$, where $e_2$ is the second standard basis element. Assume that the unit cell contains two bubbles $D_j$, $j=1,2$, each centred at $x_j$ such that

$$\begin{aligned} R_0 D_1 = D_2, \quad R_1 D_1 = D_1,\quad R_2 D_2 = D_2, \quad R_3D_1 = D_2. \end{aligned}$$

We denote the pair of bubbles by $D=D_1 \cup D_2$.

The dual lattice of $\Lambda$, denoted $\Lambda ^*$, is generated by $\alpha _1$ and $\alpha _2$ satisfying $\alpha _i\cdot l_j = 2\pi \delta _{ij}$, for $i,j = 1,2.$ Then

$$\begin{aligned} \alpha _1 = \frac{2\pi }{L}\left( \frac{1}{\sqrt{3}}, 1\right) ,~~\alpha _2 = \frac{2\pi }{L}\left( \frac{1}{\sqrt{3}}, -1 \right) . \end{aligned}$$

Fig. 2

Illustration of the dual lattice and the Brillouin zone $Y^*$

The Brillouin zone $Y^*$ is defined as the torus $Y^{*}:= {\mathbb {R}^{2}}/{\Lambda ^{*}}$ and can be represented either as

$$\begin{aligned} Y^* \simeq \left\{ s \alpha _1+ t \alpha _2 ~|~ 0 \le s,t \le 1 \right\} , \end{aligned}$$

or as the first Brillouin zone $Y_1^*$ illustrated in Fig. 2. The points

$$\begin{aligned} \alpha _1^*= \frac{2\alpha _1+\alpha _2}{3}, \quad \alpha ^*_2 = \frac{\alpha _1+2\alpha _2}{3} \end{aligned}$$

in the Brillouin zone are called Dirac points. For simplicity, in this work we will only consider the analysis around the Dirac point $\alpha ^* := \alpha ^*_1$, the main difference around $\alpha ^*_2$ is summarised in Remark 3.

We will denote the density and bulk modulus of the bubble by $\rho _b$ and $\kappa _b$, respectively. The corresponding parameters of the surrounding medium is denoted by $\rho$ and $\kappa$. Wave propagation in the bubbly honeycomb crystal is described by the following $\alpha$-quasi-periodic Helmholtz problem in Y:

$$\begin{aligned} \left\{ \begin{array} {lll} &{}\displaystyle \nabla \cdot \frac{1}{\rho } \nabla u+ \frac{\omega ^2}{\kappa } u = 0 \quad &{}\text {in} ~ Y \backslash D, \\ &{}\displaystyle \nabla \cdot \frac{1}{\rho _b} \nabla u+ \frac{\omega ^2}{\kappa _b} u = 0 \quad &{}\text {in} ~ D, \\ &{}\displaystyle u_{+} -u_{-} =0 \quad &{}\text {on}~ \partial D, \\ &{}\displaystyle \frac{1}{\rho } \frac{\partial u}{\partial \nu } \bigg |_{+} - \frac{1}{\rho _b} \frac{\partial u}{\partial \nu } \bigg |_{-} =0 \quad &{}\text {on} ~ \partial D, \\ &{}\displaystyle u(x+l)= e^{\mathrm {i}\alpha \cdot l} u(x) \quad &{} \text {for all} \ l\in \Lambda . \end{array} \right. \end{aligned}$$

(2.1)

Here, $\partial /\partial \nu$ denotes the normal derivative on $\partial D$, and the subscripts $+$ and − indicate the limits from outside and inside D, respectively. A non-trivial solution to this problem and its corresponding frequency is called a Bloch eigenfunction and a Bloch eigenfrequency. The Bloch eigenfrequencies $\omega _i^\alpha , \ i=1,2,\ldots$ with positive real part, seen as functions of $\alpha$, are known as band functions.

Let

$$\begin{aligned} v:=\sqrt{\frac{\kappa }{\rho }}, \ ~ v_b:=\sqrt{\frac{\kappa _b}{\rho _b}}, \ ~k=\frac{\omega }{v}, \ ~ k_b=\frac{\omega }{ v_b}. \end{aligned}$$

Introduce the density contrast parameter $\delta$ as

$$\begin{aligned} \delta :=\frac{\rho _b}{\rho }. \end{aligned}$$

We assume that there is a high contrast in the density while the wave speeds are comparable, i.e.,

$$\begin{aligned} \delta \ll 1 \quad \text {and} \quad v, v_b=\mathcal {O}(1). \end{aligned}$$

2.2 Quasi-Periodic Green’s Function for the Honeycomb Lattice

Define the $\alpha$-quasi-periodic Green’s function $G^{\alpha ,k}$ to satisfy

$$\begin{aligned} \Delta G^{\alpha ,k} + k^2G^{\alpha ,k} = \sum _{n \in \Lambda } \delta (x-n)e^{\mathrm {i}\alpha \cdot n}. \end{aligned}$$

Then it can be shown that $G^{\alpha ,k}$ is given by [4, 8]

$$\begin{aligned} G^{\alpha ,k}(x)= \frac{1}{|Y|}\sum _{q\in \Lambda ^*} \frac{e^{\mathrm {i}(\alpha +q)\cdot x}}{ k^2-|\alpha +q|^2}. \end{aligned}$$

(2.2)

For a given bounded domain D in Y, with Lipschitz boundary $\partial D$, the single layer potential $\mathcal {S}_D^{\alpha ,k}: L^2(\partial D) \rightarrow H_{\text {loc}}^1(\mathbb {R}^2)$ is defined by

$$\begin{aligned} {\mathcal {S}}_{D}^{\alpha ,k}[\varphi ]({x}):=\int _{\partial D} G^{\alpha ,k}( {x}-{y}) \varphi ({y}) \, \mathrm {d}\sigma ({y}),~~~{x}\in {\mathbb {R}}^{2}. \end{aligned}$$

Here, we denote by $H_{\text {loc}}^1(\mathbb {R}^2)$ the space of functions that are square integrable on every compact subset of $\mathbb {R}^2$ and have a weak first derivative that is also square integrable. The following jump relations are well-known [4, 8]:

$$\begin{aligned} \left. \frac{\partial }{\partial \nu } \mathcal {S}_{D}^{\alpha ,k}[\varphi ]\right| _{\pm }({x})=\left( \pm \frac{1}{2}I+(\mathcal {K}_{D}^{-\alpha ,k})^*\right) [\varphi ]({x}),~~~{x}\in \partial D, \end{aligned}$$

where the Neumann–Poincaré operator $(\mathcal {K}_{D}^{-\alpha ,k})^*: L^2(\partial D) \rightarrow {L}^{2}(\partial D)$ is defined as

$$\begin{aligned} (\mathcal {K}_{D}^{-\alpha ,k})^*[\varphi ](x)=\text {p.v.}\int _{\partial D}\frac{\partial }{\partial \nu _x}G^{\alpha ,k}({x}-{y}) \varphi ({y}) \, \mathrm {d}\sigma ({y}),~~~{x}\in \partial D. \end{aligned}$$

The Green’s function can be asymptotically expanded for small k as follows [4, 5]:

$$\begin{aligned} G^{\alpha ,k}= G^{\alpha ,0}+ k^{2} G_1^{\alpha } + \mathcal {O}(k^4), \qquad G_1^{\alpha }(x):=-\sum _{q\in \Lambda ^*}\frac{e^{\mathrm {i}(\alpha +q)\cdot x}}{ |\alpha +q|^{2}}, \end{aligned}$$

where the error term is uniform in $\alpha$ in a neighbourhood of $\alpha ^*$ and $x\in Y$. This leads to the following expansion of the single layer potential $\mathcal {S}_D^{\alpha ,k}: L^2(\partial D) \rightarrow H^1(\partial D)$

$$\begin{aligned} \mathcal {S}_D^{\alpha ,k}= & {} \mathcal {S}_D^{\alpha ,k} + k^{2} \mathcal {S}^{\alpha }_{D,1} + \mathcal {O}(k^4), \quad \mathcal {S}^{\alpha }_{D,1}[\phi ](x) := \int _{\partial D} G_1^{\alpha }(x-y) \nonumber \\&\times \phi (y) \, \mathrm {d}\sigma (y), \end{aligned}$$

(2.3)

where $\mathcal {O}(k^4)$ denotes an operator $L^2(\partial D) \rightarrow H^1(\partial D)$ with operator norm of order $k^4$, uniformly for $\alpha$ in a neighbourhood of $\alpha ^*$. Similarly, for the Neumann–Poincaré operator, we have

$$\begin{aligned} (\mathcal {K}_D^{-\alpha ,k})^*= & {} (\mathcal {K}_D^{-\alpha ,0})^* + k^{2} \mathcal {K}^{\alpha }_{D,1} + \mathcal {O}(k^4), \quad \mathcal {K}^{\alpha }_{D,1} [\phi ](x) := \int _{\partial D} \frac{\partial }{\partial \nu _x}G_1^{\alpha }\nonumber \\&\times (x-y) \phi (y) \, \mathrm {d}\sigma (y), \end{aligned}$$

(2.4)

where $\mathcal {O}(k^4)$ denotes an operator $L^2(\partial D) \rightarrow L^2(\partial D)$ with operator norm of order $k^4$, uniformly for $\alpha$ in a neighbourhood of $\alpha ^*$.

It is known that $\mathcal {S}_D^{\alpha ,0}: L^2(\partial D) \rightarrow H^1(\partial D)$ is invertible when $\alpha \ne 0$ [4, 8]. Let $\psi _j^{\alpha }\in L^2(\partial D)$ be given by

$$\begin{aligned} \mathcal {S}_D^{\alpha ,0}[\psi _j^{\alpha }] = \chi _{\partial D_j}\quad \hbox {on}~\partial D,\quad j=1, 2, \end{aligned}$$

(2.5)

where $\chi$ denotes the indicator function. Define the capacitance coefficient matrix $C^\alpha =(C_{ij}^\alpha )$ by

$$\begin{aligned} C_{ij}^\alpha := - \int _{\partial D_i} \psi _j^\alpha \, \mathrm {d}\sigma ,\quad i,j=1, 2. \end{aligned}$$

Using the symmetry of the honeycomb structure, it was shown in [5] that the capacitance coefficients satisfy

$$\begin{aligned} c_1^\alpha := C_{11}^\alpha = C_{22}^{\alpha }, \quad c_2^\alpha := C_{12}^{\alpha } = \overline{C_{21}^{\alpha }}, \end{aligned}$$

and

$$\begin{aligned} \nabla _\alpha c_1^\alpha \Big |_{\alpha =\alpha ^*} = 0, \quad \nabla _\alpha c_2^\alpha \Big |_{\alpha =\alpha ^*} = c \begin{pmatrix} 1\\ iu\end{pmatrix}, \end{aligned}$$

(2.6)

where we denote

$$\begin{aligned} c:=\frac{\partial c_2^{\alpha }}{\partial \alpha _1}\Big |_{\alpha =\alpha ^*}. \end{aligned}$$

In [5, Lemma 3.4] it was shown that $c\ne 0$. Note that the capacitance matrix $C^\alpha$ is written as

$$\begin{aligned} C^\alpha = \begin{pmatrix} c_1^\alpha &{}\quad c_2^\alpha \\ \overline{c_2^\alpha } &{}\quad c_1^\alpha \end{pmatrix}. \end{aligned}$$

2.3 Dirac Cone Dispersion in the Band Structure

The solution to (2.1) can be represented using the single layer potentials $\mathcal {S}_D^{\alpha , k_b}$ and $\mathcal {S}_D^{\alpha ,k}$ as follows (see, for instance, [8]):

$$\begin{aligned} u({x}) = {\left\{ \begin{array}{ll} \mathcal {S}_D^{\alpha , k_b} [\phi ]({x}),\quad {x}\in D, \\ \mathcal {S}_D^{ \alpha ,k}[\psi ]({x}),\quad {x}\in Y{\setminus }\overline{D}, \end{array}\right. } \end{aligned}$$

(2.7)

where the pair $(\phi ,\psi )\in L^2(\partial D)\times L^2(\partial D)$ satisfies

$$\begin{aligned} \mathcal {A}_\delta ^{\alpha ,\omega } \begin{pmatrix} \phi \\ \psi \end{pmatrix}=0. \end{aligned}$$

(2.8)

Here, the operator $\mathcal {A}_\delta ^{\alpha ,\omega }$ is defined by

$$\begin{aligned} \mathcal {A}_\delta ^{\alpha ,\omega } := \begin{bmatrix} \mathcal {S}^{\alpha , k_b}_D &{}\quad -\mathcal {S}^{\alpha , k}_D\\ -\frac{1}{2}I+(\mathcal {K}_{D}^{-\alpha , k_b})^* &{}\quad -\delta \left( \frac{1}{2}I+(\mathcal {K}_{D}^{-\alpha ,k})^*\right) \end{bmatrix}. \end{aligned}$$

(2.9)

It is well-known that the integral equation (2.8) has a non-trivial solution for some discrete frequencies $\omega$, and the frequencies with positive real part are the band functions $\omega _i^\alpha , \ i=1,2,\ldots$

It was shown in [5] that, for the honeycomb structure, the first two band functions $\omega _1^\alpha$ and $\omega _2^\alpha$ form a conical dispersion relation near the Dirac point $\alpha ^*$. Such a conical dispersion is referred to as a Dirac cone. More specifically, the following theorem was proved in [5] (Theorems 3.2 and 4.1 from [5]). It is worth emphasizing that the following results hold in the deep sub-wavelength regime:

Theorem 2.1

For small $\delta$, the first two band functions $\omega _j^\alpha ,j=1,2$, satisfy

$$\begin{aligned} \omega _j^\alpha = \sqrt{\frac{\delta \lambda _j^\alpha }{|D_1|}}v_b + \mathcal {O}(\delta ) \end{aligned}$$

(2.10)

uniformly for $\alpha$ in a neighbourhood of $\alpha ^*$, where $\lambda _j^\alpha , j=1,2,$ are the two eigenvalues of $C^\alpha$ and $|D_1|$ denotes the area of one of the bubbles. Moreover, for $\alpha$ close to $\alpha ^*$ and $\delta$ small enough, the first two band functions form a Dirac cone, i.e.,

$$\begin{aligned} \begin{matrix} \displaystyle \omega _1^\alpha = \omega ^*- \lambda |\alpha - \alpha ^*| \big [ 1+ \mathcal {O}(|\alpha -\alpha ^*|) \big ], \\ \displaystyle \omega _2^\alpha = \omega ^*+ \lambda |\alpha - \alpha ^*| \big [ 1+ \mathcal {O}(|\alpha -\alpha ^*|) \big ], \end{matrix} \end{aligned}$$

(2.11)

where $\omega ^*$ and $\lambda$ are independent of $\alpha$ and satisfy

$$\begin{aligned} \omega ^*= \sqrt{\frac{\delta c_1^{\alpha ^*}}{|D_1|}}v_b + \mathcal {O}(\delta ) \quad \text {and} \quad \lambda = |c|\sqrt{\delta }\lambda _0 + \mathcal {O}(\delta ), \quad \lambda _0=\frac{1}{2}\sqrt{\frac{v_b^2 }{|D_1|c_1^{\alpha ^*}}} \end{aligned}$$

as $\delta \rightarrow 0$. Moreover, the error term $\mathcal {O}(|\alpha -\alpha ^*|)$ in (2.11) is uniform in $\delta$.

In the next sections, we will investigate the asymptotic behaviour of the Bloch eigenfunctions near the Dirac points. Then we shall prove that the envelopes of the Bloch eigenfunctions satisfy a Helmholtz equation with near-zero effective refractive index and derive a two-dimensional homogenized equation of Dirac-type for the honeycomb bubbly crystal.

3 Bloch Eigenfunctions Near Dirac Points

In this section, we study the Bloch eigenfunctions in the regime close to $\alpha ^*$, following the approach of [9]. We will perform an asymptotic analysis with two small parameters. For $f,g\in C(\mathbb {R}^2,\mathbb {R})$, depending on the small parameters $\varepsilon _1,\varepsilon _2 \in \mathbb {R}$, we will use the notation

$$\begin{aligned} f(\varepsilon _1,\varepsilon _2) = \mathcal {O}(g(\varepsilon _1, \varepsilon _2)) \end{aligned}$$

to denote that there is some $K > 0$, constant in $\varepsilon _1, \varepsilon _2$, such that $|f(\varepsilon _1,\varepsilon _2)| < K|g(\varepsilon _1, \varepsilon _2)|$ for all $(\varepsilon _1, \varepsilon _2)$ in some neighbourhood of (0, 0).

3.1 Asymptotic Behaviour of the Single Layer Potential

Here we consider the single layer potential $\mathcal {S}_D^{\alpha ,k}$ near a Dirac point $\alpha = \alpha ^*$.

Lemma 3.1

Let $\psi _j^\alpha$, for $j=1,2,$ be defined by (2.5). For $j=1,2$ and for $\alpha$ near a Dirac point $\alpha ^*$, i.e. $\alpha = \alpha ^* + \varepsilon \tilde{\alpha }$ with small $\varepsilon >0$ and fixed $\tilde{\alpha }$, we have

$$\begin{aligned} \mathcal {S}_D^{\alpha ^* + \varepsilon \tilde{\alpha },k}\left[ \psi _j^{\alpha ^* + \varepsilon \tilde{\alpha }}\right] (x) = e^{\mathrm {i}\varepsilon \tilde{\alpha }\cdot x}\mathcal {S}_D^{\alpha ^*,0}\big [\xi ^{\varepsilon \tilde{\alpha }}_j\big ](x) + \mathcal {O}(\varepsilon ^2 + \varepsilon k^2), \qquad x\in \mathbb {R}, \end{aligned}$$

(3.1)

where

$$\begin{aligned} \xi ^{\varepsilon \tilde{\alpha }}_j := \left( \mathcal {S}_D^{\alpha ^*,0}\right) ^{-1}\left[ e^{-\mathrm {i}\varepsilon \tilde{\alpha }\cdot y}\chi _{\partial D_j}(y)\right] . \end{aligned}$$

Proof

As in [5], we have the Taylor expansion

$$\begin{aligned} G^{\alpha ^*+\varepsilon {\tilde{\alpha }},k}(x)&= G^{\alpha ^*,k}(x)+ \frac{1}{|Y|} \sum _{q\in \Lambda ^*} \frac{e^{\mathrm {i}(\alpha ^*+q)\cdot x}}{ k^2-|\alpha ^*+q|^2} \left( \mathrm {i}\varepsilon {\tilde{\alpha }}\cdot x +2\frac{(\alpha ^*+q)\cdot \varepsilon {\tilde{\alpha }}}{k^2-|\alpha ^*+q|^2}\right) \\&\quad +\mathcal {O}(\varepsilon ^2) = e^{\mathrm {i}\varepsilon {\tilde{\alpha }}\cdot x}\left( G^{\alpha ^*,k}(x) + G_1^{\varepsilon {\tilde{\alpha }},k}(x)\right) + \mathcal {O}(\varepsilon ^2), \end{aligned}$$

uniformly for k in a neighbourhood of 0, where

$$\begin{aligned} G_1^{\varepsilon \tilde{\alpha }, k}(x) = \frac{2}{|Y|}\sum _{q\in \Lambda ^*} \frac{e^{\mathrm {i}(\alpha ^*+q)\cdot x}(\alpha ^*+q)\cdot \varepsilon \tilde{\alpha }}{ \left( k^2-|\alpha ^*+q|^2\right) ^2}. \end{aligned}$$

We define the corresponding operator $\mathcal {S}_1^{\varepsilon \tilde{\alpha }, k}: L^2(\partial D) \rightarrow H^1(\partial D)$ as

$$\begin{aligned} \mathcal {S}_1^{\varepsilon {\tilde{\alpha }}, k}[\varphi ](x) = \int _{\partial D} G_1^{\varepsilon {\tilde{\alpha }},k}(x-y) \varphi (y) \, \mathrm {d}\sigma (y). \end{aligned}$$

Observe that, in the operator norm, $\mathcal {S}_1^{\varepsilon {\tilde{\alpha }}, k} = \mathcal {O}(\varepsilon )$, uniformly in k. With these definitions, we have the following expansion of the single layer potential for $\alpha$ close to $\alpha ^*$:

$$\begin{aligned} \mathcal {S}_D^{\alpha ^* + \varepsilon {\tilde{\alpha }},k}[\varphi ](x) = e^{\mathrm {i}\varepsilon {\tilde{\alpha }}\cdot x}\left( \mathcal {S}_D^{\alpha ^*,k} + \mathcal {S}_1^{\varepsilon {\tilde{\alpha }}, k}\right) \left[ e^{-\mathrm {i}\varepsilon {\tilde{\alpha }}\cdot y}\varphi (y)\right] (x) + \mathcal {O}(\varepsilon ^2), \end{aligned}$$

(3.2)

uniformly in k. Using the Neumann series, we can then compute

$$\begin{aligned} \psi ^{\alpha ^* + \varepsilon {\tilde{\alpha }}}_j(y)&= \left( \mathcal {S}_D^{\alpha ^*+\varepsilon {\tilde{\alpha }},0}\right) ^{-1}\left[ \chi _{\partial D_j}\right] (y) \nonumber \\&= e^{\mathrm {i}\varepsilon {\tilde{\alpha }}\cdot y}\left( I - \left( \mathcal {S}_D^{\alpha ^*,0}\right) ^{-1} \mathcal {S}_1^{\varepsilon \tilde{\alpha },0}\right) \left( \mathcal {S}_D^{\alpha ^*,0}\right) ^{-1}\nonumber \\&\quad \left[ e^{-\mathrm {i}\varepsilon {\tilde{\alpha }}\cdot x}\chi _{\partial D_j}(x)\right] (y) + \mathcal {O}(\varepsilon ^2) \nonumber \\&= e^{\mathrm {i}\varepsilon {\tilde{\alpha }}\cdot y}\left( I - \left( \mathcal {S}_D^{\alpha ^*,0}\right) ^{-1} \mathcal {S}_1^{\varepsilon {\tilde{\alpha }},0}\right) \big [\xi ^{\varepsilon {\tilde{\alpha }}}_j\big ](y) + \mathcal {O}(\varepsilon ^2),\nonumber \\&\quad \qquad y \in \partial D. \end{aligned}$$

(3.3)

Then, combining (3.2) and (3.3) and using the fact that $\mathcal {S}_D^{\alpha ,k} = \mathcal {S}_D^{\alpha ,0} +\mathcal {O}(k^2)$, uniformly for $\alpha$ in a neighbourhood of $\alpha ^*$, we obtain that

$$\begin{aligned} \mathcal {S}_D^{\alpha ^* + \varepsilon {\tilde{\alpha }},k} \left[ \psi _j^{\alpha ^* + \varepsilon {\tilde{\alpha }}} \right] (x)&= e^{\mathrm {i}\varepsilon {\tilde{\alpha }}\cdot x} \left( \mathcal {S}_D^{\alpha ^*,k} + \mathcal {S}_1^{\varepsilon {\tilde{\alpha }}, k} - \mathcal {S}_D^{\alpha ^*,k}\left( \mathcal {S}_D^{\alpha ^*,0}\right) ^{-1} \mathcal {S}_1^{\varepsilon {\tilde{\alpha }},0} \right) \\&\quad \left[ \xi ^{\varepsilon {\tilde{\alpha }}}_j\right] (x) + \mathcal {O}(\varepsilon ^2) \\&= e^{\mathrm {i}\varepsilon {\tilde{\alpha }}\cdot x}\mathcal {S}_D^{\alpha ^*,0} \big [\xi ^{\varepsilon {\tilde{\alpha }}}_j\big ](x) + \mathcal {O}(\varepsilon ^2 + \varepsilon k^2), \end{aligned}$$

which concludes the proof. $\square$

Remark 1

From the definition of $\xi ^{\varepsilon {\tilde{\alpha }}}_j$, we see that

$$\begin{aligned} \xi ^{\varepsilon {\tilde{\alpha }}}_j = \psi ^{\alpha ^*}_j + \mathcal {O}(\varepsilon ). \end{aligned}$$

(3.4)

Hence, equation (3.1) contains two terms of order $\varepsilon$: one term from the Taylor expansion of $e^{i\varepsilon {\tilde{\alpha }}\cdot x}$ and one term from the Taylor expansion of $e^{i\varepsilon {\tilde{\alpha }}\cdot y}$ inside $\xi ^{\varepsilon {\tilde{\alpha }}}_j$. Here, x and y varies on different scales: $x\in \mathbb {R}^2$ while $y\in \partial D$. Since y varies on a much smaller scale, the latter term will be negligible. This will be made precise in Section 4.

Remark 2

Equation (3.2) differs from Lemma 4.4 of [9] by the $\varepsilon$-order term $\mathcal {S}_1$ and the factor $e^{-\mathrm {i}\varepsilon {\tilde{\alpha }}\cdot y}$. This is due to small miscalculations in [9]. As seen in the final step of above proof, the term $\mathcal {S}_1$ cancels. Moreover, as in Remark 1, y varies on the small scale. In the homogenization setting defined in Section 4, the factor $e^{-\mathrm {i}\varepsilon \tilde{\alpha }\cdot y}$ only contributes to the higher orders in the expansions, and does not affect the final theorem. Therefore, despite the miscalculations in [9], the main result of [9] is still correct.

3.2 Bloch Eigenfunctions Near the Dirac Points

Here, we consider the asymptotic behaviour of the Bloch eigenfunctions near the Dirac points.

Let us assume that $\alpha$ is close to the Dirac point $\alpha ^*$, i.e., $\alpha = \alpha ^* + \varepsilon {\tilde{\alpha }}$ for small $\varepsilon >0$. Let $u^\alpha$ be the Bloch eigenfunction with the Bloch eigenfrequency $\omega ^\alpha$. In other terms, $u^\alpha$ is given by

$$\begin{aligned} u^\alpha ({x}) = {\left\{ \begin{array}{ll} \mathcal {S}_D^{\alpha , k_b} [\phi ^\alpha ]({x}),\quad {x}\in D, \\ \mathcal {S}_D^{ \alpha ,k}[\psi ^\alpha ]({x}),\quad {x}\in Y{\setminus }\overline{D}, \end{array}\right. } \end{aligned}$$

where the pair $(\phi ^\alpha ,\psi ^\alpha )\in L^2(\partial D)\times L^2(\partial D)$ satisfies

$$\begin{aligned} \mathcal {A}_\delta ^{\alpha ,\omega ^\alpha } \begin{pmatrix} \phi ^\alpha \\ \psi ^\alpha \end{pmatrix}=0. \end{aligned}$$

(3.5)

We know from Theorem 2.1 that $\omega ^\alpha = \mathcal {O}(\sqrt{\delta })$, uniformly in $\varepsilon$. So, for $\delta$ small enough and using (2.3) and (2.4), the integral equation (3.5) can be approximated by

$$\begin{aligned} {\left\{ \begin{array}{ll} \mathcal {S}_D^{\alpha ,0}[\phi ^\alpha ]- \mathcal {S}_D^{\alpha , 0}[\psi ^\alpha ] = \mathcal {O}(\delta ),\\ \left( -\frac{1}{2}I + (\mathcal {K}_D^{-\alpha ,0})^* + k_b^2 \mathcal {K}_{D,1}^{\alpha }\right) [\phi ^\alpha ]-\delta \left( \frac{1}{2}I+(\mathcal {K}_{D}^{-\alpha ,0})^*\right) [\psi ^\alpha ]= \mathcal {O}(\delta ^2). \end{array}\right. } \end{aligned}$$

(3.6)

Since $\mathcal {S}_D^{\alpha ,0}$ is invertible, and the inverse is bounded for $\alpha$ in a neighbourhood of $\alpha ^*$, the first equation in (3.6) implies

$$\begin{aligned} \psi ^\alpha = \phi ^\alpha + \mathcal {O}(\delta ), \end{aligned}$$

(3.7)

uniformly for $\alpha$ in a neighbourhood of $\alpha ^*$. Substituting (3.7) into the second equation in (3.6), we have

$$\begin{aligned} \left( -\frac{1}{2}I + (\mathcal {K}_D^{-\alpha ,0})^* + k_b^2 \mathcal {K}_{D,1}^{\alpha }\right) [\phi ^\alpha ]-\delta \left( \frac{1}{2}I+(\mathcal {K}_{D}^{-\alpha ,0})^*\right) [\phi ^\alpha ]= \mathcal {O}(\delta ^2), \end{aligned}$$

(3.8)

uniformly in $\alpha$. Since $\ker \left( -\frac{1}{2}I + (\mathcal {K}_D^{-\alpha ,0})^*\right)$ is generated by $\psi _1^\alpha$ and $\psi _2^\alpha$, which are defined by (2.5), we may write $\phi ^\alpha$ as

$$\begin{aligned} \phi ^\alpha = A \psi _1^\alpha + B\psi _2^\alpha + \mathcal {O}(\delta ), \end{aligned}$$

where $|A|+|B|=1$.

By integrating (3.8) on $\partial D_1$ and $\partial D_2$, and using the following identity [4]

$$\begin{aligned} \int _{\partial D_j} \mathcal {K}^{\alpha }_{D,1} [\phi ] \, \mathrm {d}\sigma = -\int _{D_j}\mathcal {S}_D^{\alpha ,0} [\phi ] \, \mathrm {d}x,\qquad j = 1,2, \end{aligned}$$

it follows that

$$\begin{aligned} -\frac{(\omega ^\alpha )^2 |D_1|}{v_b^2} A +\delta (A c_{1}^\alpha +B c_{2}^\alpha )= \mathcal {O}(\delta ^2),\\ -\frac{(\omega ^\alpha )^2 |D_2|}{v_b^2} B +\delta (A \overline{c_{2}^\alpha }+B c_{1}^\alpha )= \mathcal {O}(\delta ^2), \end{aligned}$$

uniformly for $\alpha$ in a neighbourhood of $\alpha ^*$. Observe that $|D_1| = |D_2|$. Thus, since we have from (2.6) that

$$\begin{aligned} c_1^{\alpha ^*+\varepsilon {\tilde{\alpha }}} = c_1^{\alpha ^*} + \mathcal {O}(\varepsilon ^2),\quad c_2^{\alpha ^*+\varepsilon {\tilde{\alpha }}} = \varepsilon c({\tilde{\alpha }}_1 - \mathrm {i}{\tilde{\alpha }}_2) + \mathcal {O}(\varepsilon ^2), \end{aligned}$$

we obtain

$$\begin{aligned} -\frac{|D_1|}{v_b^2} \left( (\omega ^\alpha )^2 - (\omega ^*)^2 \right) A + {c}\delta \varepsilon \Big ( {\tilde{\alpha }}_1 - \mathrm {i}{\tilde{\alpha }}_2\Big ) B + \mathcal {O}(\delta \varepsilon ^2)= \mathcal {O}(\delta ^2),\\ -\frac{|D_1|}{v_b^2} \left( (\omega ^\alpha )^2 - (\omega ^*)^2 \right) B + \overline{c}\delta \varepsilon \Big ( {\tilde{\alpha }}_1 + \mathrm {i}{\tilde{\alpha }}_2\Big ) A + \mathcal {O}(\delta \varepsilon ^2)= \mathcal {O}(\delta ^2). \end{aligned}$$

From (2.11) we know that

$$\begin{aligned} \omega ^\alpha + \omega ^* = 2\omega ^* + \mathcal {O}(\varepsilon \sqrt{\delta }), \quad \omega ^\alpha - \omega ^* = \mathcal {O}(\varepsilon \sqrt{\delta }), \end{aligned}$$

and from (2.10) that

$$\begin{aligned} \omega ^*= \sqrt{\frac{\delta c_1^{\alpha ^*} v_b^2}{|D_1|}} + \mathcal {O}(\delta ), \end{aligned}$$

so we arrive at

$$\begin{aligned} -2\sqrt{\frac{\delta c_1^{\alpha ^*} |D_1|}{v_b^2}} (\omega ^\alpha -\omega ^*)A + {c}\delta \varepsilon \Big ( {\tilde{\alpha }}_1 - \mathrm {i}{\tilde{\alpha }}_2\Big ) B + \mathcal {O}(\delta \varepsilon ^2)= \mathcal {O}(\delta ^2),\\ -2\sqrt{\frac{\delta c_1^{\alpha ^*} |D_1|}{v_b^2}} (\omega ^\alpha -\omega ^*)B + \overline{c}\delta \varepsilon \Big ( {\tilde{\alpha }}_1 + \mathrm {i}{\tilde{\alpha }}_2\Big ) A + \mathcal {O}(\delta \varepsilon ^2)= \mathcal {O}(\delta ^2). \end{aligned}$$

In matrix form, this reads

$$\begin{aligned}&\lambda _0 \sqrt{\delta } \begin{bmatrix} 0 &{} \varepsilon c({\tilde{\alpha }}_1 - \mathrm {i}{\tilde{\alpha }}_2)\\ \varepsilon \overline{c}( {\tilde{\alpha }}_1 + \mathrm {i}{\tilde{\alpha }}_2) &{} 0 \end{bmatrix} \begin{bmatrix} A\\ B \end{bmatrix} = (\omega ^\alpha -\omega ^*) \begin{bmatrix} A\\ B \end{bmatrix} \nonumber \\&\quad \quad + \mathcal {O}(\delta ^{1/2} \varepsilon ^2 + \delta ^{3/2}), \end{aligned}$$

(3.9)

where, as in Theorem 2.1,

$$\begin{aligned} \lambda _0 = \frac{1}{2}\sqrt{\frac{ v_b^2}{c_1^{\alpha ^*} |D_1|}}. \end{aligned}$$

Then, by solving the above eigenvalue problem, we obtain two (approximate) eigenpairs given by

$$\begin{aligned} \omega ^\alpha _\pm = \omega ^* \pm \sqrt{\delta } \lambda _0 \varepsilon |c|\cdot |{\tilde{\alpha }}| + \mathcal {O}(\delta ^{1/2} \varepsilon ^2 + \delta ^{3/2}), \end{aligned}$$

and

$$\begin{aligned} \begin{bmatrix} A_\pm \\ B_\pm \end{bmatrix} =\frac{1}{\sqrt{2}} \begin{bmatrix} \displaystyle \pm \frac{c}{|c|} \frac{{\tilde{\alpha }}_1 - \mathrm {i}{\tilde{\alpha }}_2}{|{\tilde{\alpha }}|} \\ [1.em] 1 \end{bmatrix} + \mathcal {O}(\delta + \varepsilon ^2). \end{aligned}$$

(3.10)

This implies that the Bloch eigenfunction $u_+^\alpha$ (respectively, $u_-^\alpha$) associated to the upper part $\omega _+^\alpha$ (respectively, the lower part $\omega _-^\alpha$) of the Dirac cone can be represented, for x outside D, as

$$\begin{aligned} u_\pm ^{\alpha } = A_\pm \mathcal {S}_D^{\alpha ,\omega _\pm ^\alpha /v} [\psi _1^\alpha ] + B_\pm \mathcal {S}_D^{\alpha ,\omega _\pm ^\alpha /v} [\psi _2^\alpha ] + \mathcal {O}(\delta ), \end{aligned}$$

uniformly for $\alpha$ in a neighbourhood of $\alpha ^*$. Then, by Lemma 3.1, Eq. (2.3) and the fact that $\omega _\pm ^\alpha = \mathcal {O}(\delta ^{1/2})$, uniformly for $\alpha$ in a neighbourhood of $\alpha ^*$, we have

$$\begin{aligned} u^{\alpha ^* + \varepsilon {\tilde{\alpha }}}_\pm (x)= & {} A_\pm e^{\mathrm {i}\varepsilon {\tilde{\alpha }} \cdot x} \mathcal {S}_D^{\alpha ^*,0} [\xi _1^{\varepsilon {\tilde{\alpha }}}](x) + B_\pm e^{\mathrm {i}\varepsilon {\tilde{\alpha }} \cdot x} \mathcal {S}_D^{\alpha ^*,0} [\xi _2^{\varepsilon {\tilde{\alpha }}}](x) \nonumber \\&+ \mathcal {O}(\delta + \varepsilon ^2). \end{aligned}$$

(3.11)

4 Homogenization of the Bloch Eigenfunctions Near the Dirac Points

Here, we consider the rescaled bubbly honeycomb crystal by replacing the lattice constant L with sL where $s>0$ is a small positive parameter. We then derive a homogenized equation.

We need the following lemma, which can be proved by a scaling argument:

Lemma 4.1

Let $\omega _j^\alpha ,j=1,2$, be the first two eigenvalues and $u_j^\alpha$ be the associated Bloch eigenfunctions for the bubbly honeycomb crystal with lattice constant L. Then, the bubbly honeycomb crystal with lattice constant sL has the first two Bloch eigenvalues

$$\begin{aligned} \omega _{\pm ,s}^{\alpha /s} = \frac{1}{s} \omega _\pm ^\alpha , \end{aligned}$$

and the corresponding eigenfunctions are

$$\begin{aligned} u^{\alpha /s}_{\pm ,s}(x) = u_\pm ^\alpha \left( \frac{x}{s}\right) . \end{aligned}$$

We see from the above lemma that the Dirac cone is located at the point $\alpha ^*/s$. We denote the Dirac frequency by

$$\begin{aligned} \omega _s^* = \frac{1}{s}\omega ^*. \end{aligned}$$

In the sequel, in order to simplify the presentation, we assume the following:

Assumption 4.1

The wave speed inside the bubble is equal to the one outside, i.e.,

$$\begin{aligned} v = v_b=1. \end{aligned}$$

Then the wave numbers k and $k_b$ become

$$\begin{aligned} k=\omega ^\alpha , \quad k_b=\omega ^\alpha . \end{aligned}$$

We have the following result for the Bloch eigenfunctions $u_{j,s}^{\alpha /s},j=1,2,$ for $\alpha /s$ near the Dirac points $\alpha ^*/s$.

Lemma 4.2

We have

$$\begin{aligned} u_{\pm ,s}^{\alpha ^*/s+{\tilde{\alpha }}}(x) = A_\pm e^{\mathrm {i}{\tilde{\alpha }} \cdot x} S_{1} \left( \frac{x}{s}\right) + B_\pm e^{\mathrm {i}{\tilde{\alpha }} \cdot x} S_{2} \left( \frac{x}{s}\right) + \mathcal {O}(\delta + s), \end{aligned}$$

where

$$\begin{aligned} S_j (x) = \mathcal {S}_D^{\alpha ^*,0} [\psi _j^{\alpha ^*}](x), \quad j=1,2. \end{aligned}$$

Proof

We know that $u_{\pm ,s}^{\alpha ^*/s+\tilde{\alpha }}(x) = u_\pm ^{\alpha ^* + s \tilde{\alpha }}(x/s)$, which together with (3.11) gives

$$\begin{aligned} u_{\pm ,s}^{\alpha ^*/s+{\tilde{\alpha }}}(x) = A_\pm e^{\mathrm {i}{\tilde{\alpha }} \cdot x} \mathcal {S}_D^{\alpha ^*,0} \big [\xi _1^{s{\tilde{\alpha }}}\big ]\left( \frac{x}{s}\right) + B_\pm e^{\mathrm {i}{\tilde{\alpha }} \cdot x} \mathcal {S}_D^{\alpha ^*,0} \big [\xi _2^{s{\tilde{\alpha }}}\big ]\left( \frac{x}{s}\right) + \mathcal {O}(\delta + s^2). \end{aligned}$$

By Assumption 4.1, this is valid both outside and inside D. Finally, from (3.4) we have

$$\begin{aligned} \mathcal {S}_D^{\alpha ^*,0} \big [\xi _j^{s{\tilde{\alpha }}}\big ](x) = S_j (x) + \mathcal {O}(s), \quad j=1,2, \end{aligned}$$

from which the conclusion follows. $\square$

We see that the functions $S_1$ and $S_2$ describe the microscopic behaviour of the Bloch eigenfunction $u_{\pm ,s}^{\alpha ^*/s+{\tilde{\alpha }}}$ while $A_\pm e^{\mathrm {i}{\tilde{\alpha }}\cdot x}$ and $B_\pm e^{\mathrm {i}{\tilde{\alpha }}\cdot x}$ describe the macroscopic behaviour. Now, we derive a homogenized equation near the Dirac frequency $\omega _s^*$.

Recall that the Dirac frequency of the unscaled honeycomb crystal satisfies $\omega ^* = \mathcal {O}(\sqrt{\delta })$. As in [9], in order to make the order of $\omega _s^*$ fixed when s tends to zero, we have the following assumption:

Assumption 4.2

$\delta = \mu s^2$ for some fixed $\mu >0$.

Then we have

$$\begin{aligned} \omega _s^* = \frac{1}{s} \omega ^*= \mathcal {O}(1) \quad \hbox {as } s\rightarrow 0, \end{aligned}$$

so, in what follows, we omit the subscript s in $\omega _s^*$, namely, $\omega ^*:=\omega _s^*$. Suppose the frequency $\omega$ is close to $\omega ^*$, i.e.,

$$\begin{aligned} \omega -\omega ^* = \beta \sqrt{\delta }\quad \hbox {for some constant } \beta . \end{aligned}$$

We need to find the Bloch eigenfunctions or $\tilde{\alpha }$ such that

$$\begin{aligned} \omega = \omega _{\pm ,s}^{\alpha ^*/s + {\tilde{\alpha }}}. \end{aligned}$$

We have from (3.9) and Lemmas 4.1 and 4.2 that the corresponding ${\tilde{\alpha }}$ satisfies

$$\begin{aligned} \lambda _0 \begin{bmatrix} 0 &{} c( {\tilde{\alpha }}_1 - \mathrm {i}{\tilde{\alpha }}_2)\\ \overline{c}( {\tilde{\alpha }}_1 + \mathrm {i}{\tilde{\alpha }}_2) &{} 0 \end{bmatrix} \begin{bmatrix} A_\pm \\ B_\pm \end{bmatrix} = \beta \begin{bmatrix} A_\pm \\ B_\pm \end{bmatrix} + \mathcal {O}(s). \end{aligned}$$

Thus, it is immediate to see that the macroscopic field $[{\tilde{u}}_{1}, {\tilde{u}}_{2}]^T:=[A_\pm e^{\mathrm {i}{\tilde{\alpha }}\cdot x}, B_\pm e^{\mathrm {i}{\tilde{\alpha }}\cdot x}]^T$ satisfies the system of Dirac equations as follows:

$$\begin{aligned} \lambda _0 \begin{bmatrix} 0 &{} (-c\mathrm {i})( \partial _1 - \mathrm {i}\partial _2)\\ (-\overline{c}\mathrm {i})( \partial _1 + \mathrm {i}\partial _2) &{} 0 \end{bmatrix} \begin{bmatrix} {\tilde{u}}_{1}\\ {\tilde{u}}_{2} \end{bmatrix} = \beta \begin{bmatrix} {\tilde{u}}_{1}\\ {\tilde{u}}_{2} \end{bmatrix}. \end{aligned}$$

Here, the superscript T denotes the transpose and $\partial _i$ is the partial derivative with respect to the ith variable. Note that the each component ${\tilde{u}}_j,j=1,2$, of the macroscopic field satisfies the Helmholtz equation

$$\begin{aligned} \Delta {\tilde{u}}_j + \frac{\beta ^2}{|c|^2\lambda _0^2} {\tilde{u}}_j = 0. \end{aligned}$$

(4.1)

The following is the main result on the homogenization theory for the honeycomb bubbly crystals:

Theorem 4.3

For frequencies $\omega$ close to the Dirac frequency $\omega ^*$, namely, $\omega -\omega ^* = \beta \sqrt{\delta }$, the following asymptotic behaviour of the Bloch eigenfunction $u^{\alpha ^*/s + {\tilde{\alpha }}}_s$ holds:

$$\begin{aligned} u_{s}^{\alpha ^*/s+{\tilde{\alpha }}}(x) = A e^{\mathrm {i}{\tilde{\alpha }} \cdot x} S_{1} \left( \frac{x}{s}\right) + B e^{\mathrm {i}{\tilde{\alpha }} \cdot x} S_{2} \left( \frac{x}{s}\right) + \mathcal {O}(s), \end{aligned}$$

where the macroscopic field $[{\tilde{u}}_{1}, {\tilde{u}}_{2}]^T:=[A e^{\mathrm {i}{\tilde{\alpha }}\cdot x}, B e^{\mathrm {i}{\tilde{\alpha }}\cdot x}]^T$ satisfies the two-dimensional Dirac equation

$$\begin{aligned} \lambda _0 \begin{bmatrix} 0 &{} \quad (-c\mathrm {i})( \partial _1 - \mathrm {i}\partial _2)\\ (-\overline{c}\mathrm {i})( \partial _1 + \mathrm {i}\partial _2) &{} \quad 0 \end{bmatrix} \begin{bmatrix} {\tilde{u}}_{1}\\ {\tilde{u}}_{2} \end{bmatrix} = \frac{\omega -\omega ^*}{\sqrt{\delta }} \begin{bmatrix} {\tilde{u}}_{1}\\ {\tilde{u}}_{2} \end{bmatrix}, \end{aligned}$$

which can be considered as a homogenized equation for the honeycomb bubbly crystal while the microscopic fields $S_1$ and $S_2$ vary on the scale of s.

Remark 3

Theorem 4.3 is valid around the Dirac point $\alpha ^* = \alpha _1^*$. Around the other Dirac point, analogous arguments show that Theorem 4.3 is valid with all quantities instead defined using $\alpha ^*=\alpha _2^*$ and the macroscopic field now satisfying

$$\begin{aligned} \lambda _0 \begin{bmatrix} 0 &{} \quad (-c\mathrm {i})( \partial _1 + \mathrm {i}\partial _2)\\ (-\overline{c}\mathrm {i})( \partial _1 - \mathrm {i}\partial _2) &{}\quad 0 \end{bmatrix} \begin{bmatrix} {\tilde{u}}_{1}\\ {\tilde{u}}_{2} \end{bmatrix} = \frac{\omega -\omega ^*}{\sqrt{\delta }} \begin{bmatrix} {\tilde{u}}_{1}\\ {\tilde{u}}_{2} \end{bmatrix}. \end{aligned}$$

5 Numerical Illustrations

In this section, we illustrate the main result of this paper, namely Theorem 4.3, in the case of circular bubbles. We do this by numerically computing the eigenmodes close to the Dirac points for the honeycomb lattice. For comparison, in Section 5.2, we also compute the eigenmodes in the case of a square lattice of bubbles. This will show the necessity of having a honeycomb lattice, instead of the simpler square lattice, in order to achieve near-zero index in the sub-wavelength regime. This also serves to illustrate and numerically verify the conclusions from [9].

The eigenmodes are computed by discretising the operator $\mathcal {A}_\delta ^{\alpha ,\omega }$ from Eq. (2.9) using the multipole method as described in [5, 6]. All computations are made for circular bubbles with radius $R= 0.2$. Moreover, the material parameters are $\rho = \kappa = 1000$, $\rho _b = \kappa _b = 1$, which gives $\delta = 10^{-3}$, $v=1$, and $v_b=1$.

For simplicity, we will perform the computations with the scaling $s=1$. Observe that Theorem 4.3 is valid for small s, which will instead be achieved using the rescaling properties of the eigenmodes. From the proof of Lemma 4.2 we have $u_{\pm ,s}^{\alpha ^*/s+{\tilde{\alpha }}}(x) = u_\pm ^{\alpha ^* + s {\tilde{\alpha }}}(x/s)$. Therefore, a small scaling factor s corresponds to choosing quasi-periodicities $\alpha$ sufficiently close to $\alpha ^*$ and choosing a large range of x. This justifies the choice $s=1$, and allows us to study the limiting behaviour without changing the geometry of the differential equation.

5.1 Honeycomb Lattice

Recall from (2.7) that the eigenmodes can be expressed as

$$\begin{aligned} u({x}) = {\left\{ \begin{array}{ll} \mathcal {S}_D^{ \alpha ,k}[\psi ]({x}),\quad {x}\in Y{\setminus }\bar{D},\\ \mathcal {S}_D^{\alpha , k_b} [\phi ]({x}),\quad {x}\in D, \end{array}\right. } \end{aligned}$$

(5.1)

where $\mathcal {A}_\delta ^{\alpha ,\omega } \left( {\begin{matrix} \phi \\ \psi \end{matrix}}\right) =0.$ Then $(\phi ,\psi )$ can be numerically computed as an eigenvector of the discretised operator $\mathcal {A}_\delta ^{\alpha ,\omega }$ corresponding to the eigenvalue 0. Moreover, u(x) can be computed by (5.1), extended quasi-periodically to the whole space. We will consider the eigenmodes at frequencies $\omega$ with $\omega -\omega ^* = \beta \sqrt{\delta }$.

First, we consider the small-scale behaviour of the eigenmodes. The small scale corresponds to $e^{\mathrm {i}{\tilde{\alpha }}\cdot x} \approx 1$, so Theorem 4.3 shows that the eigenfunctions are given by

$$\begin{aligned} u_{\pm }^{\alpha ^*+{\tilde{\alpha }}}(x) = A S_{1} (x) + B S_{2} (x) + \mathcal {O}(s). \end{aligned}$$

Equation (3.10) shows that $A = A_\pm = \pm \frac{1}{\sqrt{2}}e^{\mathrm {i}(\theta _c+\theta )}$ and $B = B_\pm = \frac{1}{\sqrt{2}}$, where $\theta _c$ and $\theta$ are the arguments of c and ${\tilde{\alpha }}$, respectively, and the sign coincides with the sign of $\beta$. To pick a unique eigenmode, we choose $\theta =0$, i.e., the quasi-periodicity ${\tilde{\alpha }} = \left( {\begin{matrix} {\tilde{\alpha }}_1 \\ 0 \end{matrix}}\right)$.

Figure 3 shows the function

$$\begin{aligned} u^*(x) = A_- S_{1} (x) + B_- S_{2} (x), \end{aligned}$$

which is the first Dirac eigenmode in the limit $\beta \rightarrow 0^-$. It can be seen that the eigenmode is highly oscillating and oscillates between $-1$ and 1 within one hexagon of bubbles. Moreover, this eigenmode has no large-scale oscillations. The second eigenmode, corresponding to $\beta \rightarrow 0^+$, has the same qualitative features.

Next, we consider the large-scale behaviour of the eigenmodes. Figure 4 shows the real part of the first eigenmode $u_{-}^{\alpha ^*+{\tilde{\alpha }}}$ for $\beta = 8\cdot 10^{-3}$. It can be seen that the eigenmode is oscillating with a low frequency in the large scale. Moreover, it is clear that the eigenmode consists of two superimposed fields, corresponding to the parts $Ae^{\mathrm {i}{\tilde{\alpha }}\cdot x} S_{1} (x)$ and $B e^{\mathrm {i}{\tilde{\alpha }}\cdot x}S_{2} (x)$. These fields are phase-shifted, due to the factor $e^{\mathrm {i}(\theta _c + \theta )}$ in A.

To demonstrate the near-zero effective refractive index property of the bubbly honeycomb crystal, we consider the large-scale spatial oscillation frequency close to the Dirac points. From (4.1), we know that for frequencies $\omega = \omega ^*+ \varepsilon$ close to the Dirac frequency, each component $\tilde{u}_j, j=1,2,$ of the macroscopic field oscillates at a spatial frequency

$$\begin{aligned} f = \frac{|\varepsilon |}{|c|\lambda _0\sqrt{\delta }}. \end{aligned}$$

To verify this relation, we compute the large-scale spatial frequency of the eigenmodes for $\varepsilon$ in the range $\varepsilon \in [-\,0.01,0.01]$, as shown in Fig. 5. It can be seen that the relation is linear for $\varepsilon$ close to 0. This verifies (4.1), and shows that close to the Dirac frequency the macroscopic behaviour of the honeycomb crystal can be described as a near-zero refractive index material. However, we emphasize that the homogenization approach is valid only in the large scale; this will fail to capture the small-scale oscillations. Also, we emphasize the counter-intuitive result that despite not being located at $\Gamma$, the eigenmodes close to the Dirac point show close to zero phase change across the crystal. Indeed, this is a consequence of the near-zero property and the zero phase change of the microscopic field across a hexagon.

Fig. 3

Small-scale behaviour of the first Bloch eigenfunction $u^*$

Fig. 4

Real part of first Bloch eigenfunction of the honeycomb lattice shown over many unit cells

Fig. 5

Spatial frequency of the macroscopic functions ${\tilde{ u}}_j$ for different frequency shifts $\varepsilon :=\beta \sqrt{\delta }$

5.2 Square Lattice

In this section, we perform the same numerical experiments as in Section 5.1 but for the case of a square lattice of bubbles. We begin by recalling the main result from [9]. Consider now a square lattice with unit cell Y as depicted in Fig. 6. The corresponding dispersion relation has a bandgap between the first and second bands, and the critical frequency $\omega ^*$ of the first band is achieved at the symmetry point $\alpha ^* = M = (\pi ,\pi )$ in the Brillouin zone.

Fig. 6

Illustration of the square lattice crystal and quantities in the fundamental domain Y

Now let $S_D^{\alpha ,k}$ denote the single-layer potential for the square lattice, defined analogously as in Section 2.2 but with $\Lambda$ and $\Lambda ^*$ being the square lattice and dual square lattice, respectively. Define the function

$$\begin{aligned} S(x) = \mathcal {S}^{\alpha ,0}_D\left[ \psi ^{\alpha ^*}\right] (x), \quad x\in \mathbb {R}^2, \end{aligned}$$

where $\psi ^{\alpha } = \left( \mathcal {S}^{\alpha ,0}_D\right) ^{-1}\left[ \chi _{\partial D}\right]$.

We now consider the rescaled crystal with unit cell $sY, s>0$. Again, we assume that the order of $\omega ^*$ is fixed, i.e., $\delta = \mu s^2,$ for some $\mu >0$. Then the eigenmodes of the rescaled square crystal are given in the following theorem, which is the analogue of Theorem 4.3 for the case of a square lattice:

Theorem 5.1

[9] For frequencies $\omega$ close to the critical frequency $\omega ^*$, namely, $\left( \omega ^*\right) ^2 - \omega ^2 = \mathcal {O}(s^2)$, the following asymptotic behaviour of the Bloch eigenfunction $u^{\alpha ^*/s + {\tilde{\alpha }}}_s$ holds:

$$\begin{aligned} u_{s}^{\alpha ^*/s+{\tilde{\alpha }}}(x) = e^{\mathrm {i}{\tilde{\alpha }} \cdot x} S\left( \frac{x}{s}\right) + \mathcal {O}(s), \end{aligned}$$

where the macroscopic field ${\tilde{u}}:= e^{\mathrm {i}{\tilde{\alpha }}\cdot x}$ satisfies the Helmholtz equation

$$\begin{aligned} \Delta {\tilde{u}} + \frac{(\omega ^*)^2-\omega ^2}{\delta {\tilde{\lambda ^2}}} {\tilde{u}} = 0, \end{aligned}$$

(5.2)

which can be considered as a homogenized equation for the square bubbly crystal, while the microscopic field S varies on the scale of s.

The isotropic form of the macroscopic Eq. (5.2) follows since D is a circle, and an expression for ${\tilde{\lambda }}$ is given in [9].

Fig. 7

Real part of small-scale behaviour of the Bloch eigenfunction of the square lattice at $\alpha =\alpha ^*$ (imaginary part close to 0)

Fig. 8

Real part of Bloch eigenfunction of the square lattice shown over many unit cells

Fig. 9

Spatial frequency of the envelope function ${\tilde{u}}$ for different frequency shifts $\varepsilon$ in the case of a square lattice

We now compute the eigenmodes of the square crystal close to the critical frequency, namely, $\omega =\omega ^*-\varepsilon$. The small-scale behaviour of the eigenmodes, i.e., the function $S\left( \frac{x}{s}\right)$, is shown in Fig. 7. It can be seen that the function oscillates at the scale of the bubbles. Next, the large-scale behaviour is considered. Figure 8 shows the Bloch eigenfunction over many unit cells for ${\tilde{\alpha }} = \left( {\begin{matrix} {\tilde{\alpha }}_1 \\ 0 \end{matrix}}\right)$ and $\omega ^*-\omega = 6\cdot 10^{-5}$. Similarly as in the case of a honeycomb crystal, the eigenfunction varies at the large scale with a low-frequency macroscopic field ${\tilde{u}}$.

To illustrate the macroscopic Eq. (5.2), the spatial frequency f of the macroscopic field is computed for $\varepsilon$ in the range $\varepsilon \in [0,0.01]$. Observe that due to the bandgap above $\omega ^*$, no eigenmodes exists for $\varepsilon < 0$. Figure 9 shows the spatial frequency f for different $\varepsilon$. This figure shows that f scales like $\sqrt{\varepsilon }$ for small $\varepsilon > 0$. This is consistent with (5.2), from which we expect a spatial frequency

$$\begin{aligned} f = \sqrt{\frac{(\omega ^*)^2-\omega ^2}{{\tilde{\lambda }}^2\delta }} = \sqrt{\varepsilon }\left( \sqrt{\frac{\omega ^*+\omega }{{\tilde{\lambda }}^2\delta }}\right) . \end{aligned}$$

In summary, Figs. 4 and 8 show that both in the cases of a honeycomb lattice and a square lattice, the eigenmodes have a periodic small-scale oscillation and a large-scale macroscopic oscillation. However, Figs. 5 and 9 show that the spatial frequency of the macroscopic fields have different asymptotic behaviours close to the critical frequency, due to the fact that the square lattice cannot be mapped to a zero-index effective material.

6 Concluding Remarks

In this paper we have derived, for the first time, the equation governing wave propagation in a honeycomb crystal of sub-wavelength resonators near the Dirac points. We have decomposed the Bloch eigenfunctions as the sum of two eigenmodes. The effective equation for the envelope of each of these eigenmodes is of Helmholtz-type with near-zero refractive index. Furthermore, we have shown that the two envelopes are phase-shifted and satisfy a system of Dirac equations. A comparison with a square lattice crystal shows the great potential of using honeycomb crystals of sub-wavelength resonators as near-zero materials. In a forthcoming work, we plan to study topological phenomena in bubbly time-dependent crystals and derive their effective properties. These materials may exhibit a Dirac cone near the point $\Gamma$ and therefore, may exhibit finite acoustic impedance and a high transmittance [14]. We also plan to mathematically analyse wave propagation phenomena in near-zero refractive index materials.