A high-frequency homogenization approach near the Dirac points in bubbly honeycomb crystals

TL;DR

This paper develops a high-frequency homogenization method to analyze Bloch eigenfunctions near Dirac points in bubbly honeycomb crystals, revealing a near-zero effective refractive index and detailed eigenmode decomposition.

Contribution

It introduces a novel homogenization approach that decomposes eigenmodes into slowly varying and oscillatory parts, demonstrating a near-zero refractive index near Dirac points in sub-wavelength metamaterials.

Findings

Eigenfunctions decompose into two components: slowly varying and highly oscillating.

Slowly varying components satisfy Dirac equations.

Near-zero effective refractive index is achieved at Dirac points.

Abstract

In [H. Ammari et al., Honeycomb-lattice Minnaert bubbles. arXiv:1811.03905], the existence of a Dirac dispersion cone in a bubbly honeycomb phononic crystal 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 proves 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.

Equations

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

$$Y:=\left\{s{l}_1+t{l}_2|0\le s,t\le 1\right\}.$$

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

$$R_{1}x=Rx+l_1,R_{2}x=Rx+2l_1,R_{0}x=2x_0-x,$$

$$R_0{D}_1=D_2,R_1{D}_1=D_1,R_2{D}_2=D_2,R_3{D}_1=D_2.$$

$${\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).$$

$$Y^{\ast }\simeq \left\{s{\alpha }_1+t{\alpha }_2|0\le s,t\le 1\right\},$$

$${\alpha }_1^{\ast }=\frac{2{\alpha }_1+{\alpha }_2}{3},{\alpha }_2^{\ast }=\frac{{\alpha }_1+2{\alpha }_2}{3},$$

$$\left\{\begin{array}[]{lll}&\displaystyle\nabla\cdot\frac{1}{\rho}\nabla u+\frac{\omega^{2}}{\kappa}u=0&\text{in}\ Y\backslash D,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&\displaystyle\nabla\cdot\frac{1}{\rho_{b}}\nabla u+\frac{\omega^{2}}{\kappa_{b}}u=0&\text{in}\ D,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&\displaystyle u_{+}-u_{-}=0&\text{on}\ \partial D,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&\displaystyle\frac{1}{\rho}\frac{\partial u}{\partial\nu}\bigg{|}_{+}-\frac{1}{\rho_{b}}\frac{\partial u}{\partial\nu}\bigg{|}_{-}=0&\text{on}\ \partial D,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&u(x+l)=e^{\mathrm{i}\mkern 1.0mu\alpha\cdot l}u(x)&\text{for all}\ l\in\Lambda.\end{array}\right.$$

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

$$\delta :=\frac{{\rho }_b}{\rho }.$$

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

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

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

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

$${\frac{\partial }{\partial \nu }{\mathcal{S}}_D^{\alpha ,k}[\phi ]|}_{\pm }(x)=\left(\pm \frac{1}{2}I+({\mathcal{K}}_D^{-\alpha ,k}{)}^{\ast }\right)[\phi ](x),x\in \partial D,$$

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

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

$${\mathcal{S}}_D^{\alpha ,k}={\mathcal{S}}_D^{\alpha ,k}+k^2{\mathcal{S}}_{D,1}^{\alpha }+\mathcal{O}(k^4),{\mathcal{S}}_{D,1}^{\alpha }[\varphi ](x):={\int }_{\partial D}{G}_1^{\alpha }(x-y)\varphi (y)\mathrm{d}\sigma (y),$$

$$({\mathcal{K}}_D^{-\alpha ,k}{)}^{\ast }=({\mathcal{K}}_D^{-\alpha ,0}{)}^{\ast }+k^2{\mathcal{K}}_{D,1}^{\alpha }+\mathcal{O}(k^4),{\mathcal{K}}_{D,1}^{\alpha }[\varphi ](x):={\int }_{\partial D}\frac{\partial }{\partial {\nu }_x}{G}_1^{\alpha }(x-y)\varphi (y)\mathrm{d}\sigma (y),$$

$${\mathcal{S}}_D^{\alpha ,0}[{\psi }_j^{\alpha }]={\chi }_{\partial D_j}\text{mbox}on\partial D,j=1,2,$$

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

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

$$\nabla_{\alpha}c_{1}^{\alpha}\Big{|}_{\alpha=\alpha^{*}}=0,\quad\nabla_{\alpha}c_{2}^{\alpha}\Big{|}_{\alpha=\alpha^{*}}=c\begin{pmatrix}1\\ -\mathrm{i}\mkern 1.0mu\end{pmatrix},$$

$$c:=\frac{\partial c_{2}^{\alpha}}{\partial\alpha_{1}}\Big{|}_{\alpha=\alpha^{*}}.$$

$$C^{\alpha }=\left(\begin{array}{cc}{c}_1^{\alpha }& c_2^{\alpha }\\ \overline{c_2^{\alpha }}& c_1^{\alpha }\end{array}\right).$$

$$u(x)=\{\begin{array}{l}{\mathcal{S}}_D^{\alpha ,k_b}[\varphi ](x),x\in D,\\ {\mathcal{S}}_D^{\alpha ,k}[\psi ](x),x\in Y\setminus \overline{D},\end{array}$$

$${\mathcal{A}}_{\delta }^{\alpha ,\omega }\left(\begin{array}{c}\varphi \\ \psi \end{array}\right)=0.$$

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

$${\omega }_j^{\alpha }=\sqrt{\frac{\delta {\lambda }_j^{\alpha }}{|D_1|}}{v}_b+\mathcal{O}(\delta ),$$

Full text

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

Keywords.  Honeycomb lattice, Dirac cone, Dirac equation, bubble, Minnaert resonance, sub-wavelength bandgap, near-zero effective index.

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, 22, 31, 40].

The first near-zero refractive index phononic crystal was theoretically demonstrated in [31], where the effective mass density and reciprocal bulk modulus were shown to vanish simultaneously. This near-zero effective refractive index is a consequence of the existence of a Dirac dispersion cone in the dispersion relation of the material. 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, 20, 21]. At high frequencies, the Bloch eigenmodes will oscillate on the microscale of the metamaterial, suggesting that a homogeneous description of the material is overly simplified. Nevertheless, as will be shown in this paper, an effective equation for the envelopes of these Bloch eigenmodes can be derived.

Metamaterials with Dirac singularities have been experimentally and numerically studied in [35, 38, 39]. 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, 17, 26, 34, 36]. 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, 33]. 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 [23, 25, 28, 29, 30, 32, 41] have been rigorously explained in [3, 6, 7, 5, 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. This generalizes the result obtained in [19], for the Schrödinger equation, to wave propagation in sub-wavelength resonant structures. The main result 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, and demonstrate the near-zero effective refractive index of the honeycomb 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 Figure 1. Define the lattice $\Lambda$ generated by the lattice vectors

[TABLE]

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

[TABLE]

Define the three points $x_{0},x_{1},$ and $x_{2}$ as

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

The Brillouin zone $Y^{*}$ is defined as the torus Y^{*}:={\raisebox{1.99997pt}{\mathbb{R}^{2}}\left/\raisebox{-1.99997pt}{\Lambda^{*}}\right.} and can be represented either as

[TABLE]

or as the first Brillouin zone $Y_{1}^{*}$ illustrated in Figure 2. The points

[TABLE]

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$:

[TABLE]

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,...$ with positive real part, seen as functions of $\alpha$, are known as band functions.

Let

[TABLE]

Introduce the density contrast parameter $\delta$ as

[TABLE]

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

[TABLE]

2.2 Quasi-periodic Green’s function for the honeycomb lattice

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

[TABLE]

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

[TABLE]

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_{\textrm{loc}}^{1}(\mathbb{R}^{2})$ is defined by

[TABLE]

Here, we denote by $H_{\textrm{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]:

[TABLE]

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

[TABLE]

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

[TABLE]

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)$

[TABLE]

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

[TABLE]

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\neq 0$ [4, 8]. Let $\psi_{j}^{\alpha}\in L^{2}(\partial D)$ be given by

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

where we denote

[TABLE]

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

[TABLE]

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]):

[TABLE]

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

[TABLE]

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

[TABLE]

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,...$

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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 $\epsilon_{1},\epsilon_{2}\in\mathbb{R}$, we will use the notation

[TABLE]

to denote that there is some $K>0$, constant in $\epsilon_{1},\epsilon_{2}$, such that $|f(\epsilon_{1},\epsilon_{2})|<K|g(\epsilon_{1},\epsilon_{2})|$ for all $(\epsilon_{1},\epsilon_{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^{*}+\epsilon\tilde{\alpha}$ with small $\epsilon>0$ and fixed $\tilde{\alpha}$, we have

[TABLE]

where

[TABLE]

Proof.

As in [5], we have the Taylor expansion

[TABLE]

uniformly for $k$ in a neighbourhood of [math], where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

which concludes the proof. ∎

Remark 1**.**

From the definition of $\xi^{\epsilon\tilde{\alpha}}_{j}$, we see that

[TABLE]

Hence, equation (3.1) contains two terms of order $\epsilon$: one term from the Taylor expansion of $e^{i\epsilon\tilde{\alpha}\cdot x}$ and one term from the Taylor expansion of $e^{i\epsilon\tilde{\alpha}\cdot y}$ inside $\xi^{\epsilon\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 [9, Lemma 4.4] by the $\epsilon$-order term $\mathcal{S}_{1}$ and the factor $e^{-\mathrm{i}\mkern 1.0mu\epsilon\tilde{\alpha}\cdot y}$, which, in general, do not cancel. This is due to small miscalculations in [9]. Nevertheless, under the homogenization setting defined in Section 4 and [9, Section 5], these factors only contributes to the higher orders in the expansions, and does not affect the final theorems. Therefore, despite the miscalculations in [9], the main result, [9, Theorem 5.4] 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^{*}+\epsilon\tilde{\alpha}$ for small $\epsilon>0$. Let $u^{\alpha}$ be the Bloch eigenfunction with the Bloch eigenfrequency $\omega^{\alpha}$. In other terms, $u^{\alpha}$ is given by

[TABLE]

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

[TABLE]

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

[TABLE]

uniformly for $\alpha$ in a neighbourhood of $\alpha^{*}$. 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

[TABLE]

uniformly in $\alpha$. Substituting (3.7) into the second equation in (3.6), we have

[TABLE]

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

[TABLE]

uniformly in $\alpha$, where $|A|+|B|=1$.

By integrating (3.8) on $\partial D_{1}$ and $\partial D_{2}$, and using the following identity [4]

[TABLE]

it follows that

[TABLE]

uniformly for $\alpha$ in a neighbourhood of $\alpha^{*}$. Observe that $|D_{1}|=|D_{2}|$. Thus, since we have from (2.6) that

[TABLE]

we obtain

[TABLE]

From (2.11) we know that

[TABLE]

and from (2.10) that

[TABLE]

so we arrive at

[TABLE]

In matrix form, this reads

[TABLE]

where, as in Theorem 2.1,

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and the corresponding eigenfunctions are

[TABLE]

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

[TABLE]

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.,

[TABLE]

Then the wave numbers $k$ and $k_{b}$ become

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

from which the conclusion follows. ∎

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}\mkern 1.0mu\tilde{\alpha}\cdot x}$ and $B_{\pm}e^{\mathrm{i}\mkern 1.0mu\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 assume that

Assumption 4.2**.**

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

Then we have

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

where the macroscopic field $[\tilde{u}_{1},\tilde{u}_{2}]^{T}:=[Ae^{\mathrm{i}\mkern 1.0mu\tilde{\alpha}\cdot x},Be^{\mathrm{i}\mkern 1.0mu\tilde{\alpha}\cdot x}]^{T}$ satisfies the two-dimensional Dirac equation

[TABLE]

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 $\alpha^{*}=\alpha_{2}^{*}$ and the macroscopic field now satisfying

[TABLE]

where $c,S_{1}$ and $S_{2}$ are now defined using $\alpha^{*}=\alpha_{2}^{*}$.

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 equation (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

[TABLE]

where $\mathcal{A}_{\delta}^{\alpha,\omega}\left(\begin{smallmatrix}\phi\\ \psi\end{smallmatrix}\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 [math]. 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}\mkern 1.0mu\tilde{\alpha}\cdot x}\approx 1$, so Theorem 4.3 shows that the eigenfunctions are given by

[TABLE]

Equation (3.10) shows that $A=A_{\pm}=\pm\frac{1}{\sqrt{2}}e^{\mathrm{i}\mkern 1.0mu(\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{smallmatrix}\tilde{\alpha}_{1}\\ 0\end{smallmatrix}\right)$.

Figure 3 shows the function

[TABLE]

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}\mkern 1.0mu\tilde{\alpha}\cdot x}S_{1}(x)$ and $Be^{\mathrm{i}\mkern 1.0mu\tilde{\alpha}\cdot x}S_{2}(x)$. These fields are phase-shifted, due to the factor $e^{\mathrm{i}\mkern 1.0mu(\theta_{c}+\theta)}$ in $A$.

To demonstrate the near-zero effective refractive index 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^{*}+\epsilon$ close to the Dirac frequency, each component $\tilde{u}_{j},j=1,2,$ of the macroscopic field oscillates at a spatial frequency

[TABLE]

To verify this relation, we compute the large-scale spatial frequency of the eigenmodes for $\epsilon$ in the range $\epsilon\in[-0.01,0.01]$, as shown in Figure 5. It can be seen that the relation is linear for $\epsilon$ close to [math]. 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.

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 Figure 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.

Let now $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

[TABLE]

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:

[TABLE]

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

[TABLE]

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 equation (5.2) follows since $D$ is a circle, and an expression for $\tilde{\lambda}$ is given in [9].

We now compute the eigenmodes of the square crystal close to the critical frequency, namely, $\omega=\omega^{*}-\epsilon$. The small-scale behaviour of the eigenmodes, i.e., the function $S\left(\frac{x}{s}\right)$, is shown in Figure 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{smallmatrix}\tilde{\alpha}_{1}\\ 0\end{smallmatrix}\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 equation (5.2), the spatial frequency $f$ of the macroscopic field is computed for $\epsilon$ in the range $\epsilon\in[0,0.01]$. Observe that due to the bandgap above $\omega^{*}$, no eigenmodes exists for $\epsilon<0$. Figure 9 shows the spatial frequency $f$ for different $\epsilon$. This figure shows that $f$ scales like $\sqrt{\epsilon}$ for small $\epsilon>0$. This is consistent with (5.2), from which we expect a spatial frequency

[TABLE]

In summary, Figures 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, Figures 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.