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

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 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 Dirac equation. Our results in this paper demonstrate 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 subwavelength 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.


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 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 Γ 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].
In [9], the existence of a Dirac dispersion cone near the symmetry point K in a bubbly honeycomb phononic crystal comprised of bubbles of arbitrary shape is shown. In this paper, we demonstrate near the point K the near-zero effective property of a bubbly honeycomb crystal at the deep sub-wavelength scale. In order to achieve high transmittance, bubbly crystals with time-dependent material parameters should be used [14] to turn the point K to the centre Γ of the Brillouin zone.
Metamaterials with Dirac singularities have been experimentally and numerically studied in [38,39,35]. 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].
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 subwavelength 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 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 away from the centre Γ 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 envelope 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 in the vicinity of M 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 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 Dirac equation 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 [19] for the Schrödinger equation. 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. 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

Problem formulation
In this section, we describe the honeycomb lattice and state the main results of [5]. We consider a two-dimensional infinite honeycomb crystal in two dimensions depicted in 1. Define a hexagonal lattice Λ with lattice vectors: Denote by Y a fundamental domain of the given lattice. Here, we take Define the three points x 0 , x 1 and x 2 as We will assume that each bubble in the crystal has a three-fold rotational symmetry and that each pair of adjacent bubbles has a two-fold rotational symmetry. More precisely, let R 1 and R 2 be the rotations by − 2π 3 around x 1 and x 2 , respectively, and let R 0 be the rotation around x 0 by π. These rotations can be written as where R is the rotation by − 2π 3 around the origin. Assume that the unit cell contains two bubbles D j , j = 1, 2, each centred at x j such that Denote the bubble dimer by D = D 1 ∪ D 2 .
x 0 x 2 Figure 1: Illustration of the bubbly honeycomb crystal and quantities in the fundamental domain Y . Figure 2: Illustration of the dual lattice and the Brillouin zone Y * .
The dual lattice of Λ, denoted Λ * , is generated by α 1 and α 2 satisfying α i · l j = 2πδ ij for i, j = 1, 2. Then The Brillouin zone Y * is defined as the torus Y * := R 2 / Λ * and can be represented either as or as the first Brillouin zone Y * 1 . The points in the Brillouin zone are called Dirac points.
Wave propagation in the bubbly honeycomb crystal is described by the following α-quasi-periodic Helmholtz problem in Y : (2.1) Here, ∂ ∂ν denotes the normal derivative on ∂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 the Bloch eigenfunction and the Bloch eigenfrequency. Let Introduce the density contrast parameter δ as We assume that there is a high contrast in the density while the wave speeds are comparable, i.e., δ ≪ 1 and v, v b = O(1).

Quasi-periodic Green's function for the honeycomb lattice
Define the α−quasi-periodic Green's function G α,k to satisfy Then it can be shown that G α,k is given by [4,8] For a given bounded domain D in Y , with Lipschitz boundary ∂D, the single layer potential of the density function ϕ ∈ L 2 (∂D) is defined by The following jump relations are well-known [4,8]: It is known that S α,0 D : where χ denotes the indicator function. Define the capacitance coefficient matrix Using the symmetry of the honeycomb structure, it was shown in [5] that the capacitance coefficients satisfy where we denote c = ∂c α 2 ∂α 1 α=α * . Note that the capacitance matrix C α is written as (2.8)

Dirac cone dispersion in the band structure
The solution to (2.1) can be represented using the single layer potentials S α,k b D and S α,k D as follows (see for example [8]): Here, the operator A α,ω δ is defined by It is well-known that the integral equation (2.10) has a non-trivial solution for some discrete frequencies ω. These can be viewed as the characteristic values of the operatorvalued analytic function A α,ω δ (with respect to ω); see [8]. It can be shown that ω = 0 is a characteristic value of A α,ω has non-trivial kernel of two dimensions, which is generated by where ψ α 1 and ψ α 2 are as in (2.6). Then Gohberg-Sigal theory [8] tells us that there exists characteristic values ω α j = ω α j (δ), j = 1, 2 of A α,ω δ such that ω α j (0) = 0 and ω α j depends on δ continuously.
It was shown in [5] that, for the honeycomb structure, the first two characteristic values ω α 1 and ω α 2 form a conical dispersion relation near the Dirac point α * . Such a conical dispersion is referred to as a Dirac cone. More specifically, the following theorem was proved in [5]. It worth emphasizing that the following results hold in the deep sub-wavelength regime.
In the next sections, we will investigate the asymptotic behaviour of the Bloch eigenfunctions near Dirac points. Then we shall derive a homogenized equation.

Bloch eigenfunctions near Dirac points
In this section, we study the Bloch eigenfunctions in the regime close to α * , following the approach of [9].

Asymptotic behaviour of the Green's function
In this section, we consider G α,k near a Dirac point α = α * .

3
. It is easily seen that We will use the symmetry of the problem to show that the absolutely convergent series vanishes, which proves the lemma. We define three sub-lattices Λ * 1 , Λ * 2 , and Λ * 3 in the following way. Let Λ * 1 be generated by the vectors q 1 = 2α 1 + α 2 and q 2 = α 1 + 2α 2 , and define the other two sub-lattices as It can be shown that Λ * 1 , Λ * 2 , and Λ * 3 is a disjoint partition of Λ * and that RΛ * 1 = Λ * 1 . From these facts, we find Therefore, form the fact that I + R + R 2 = 0, it follows that which proves the claim.

Bloch eigenfunctions near the Dirac points
Here we consider the asymptotic behaviour of the Bloch eigenfunctions near the Dirac points. Let us assume that α is close to the Dirac point α * , i.e., α = α * + ǫα for small ǫ > 0. Let u α be the Bloch eigenfunction with the Bloch eigenfrequency ω α . In other terms, u α is given by where the pair (φ α , ψ α ) ∈ L 2 (∂D) × L 2 (∂D) satisfies In the sequel, in order to simplify the presentation, we assume the following: Then the wave numbers k and k b become We know from Theorem 2.1 that ω α = O( √ δ). So, for δ small enough, the integral equation (3.5) can be approximated by Then, since S α,0 D is invertible when α = 0, the first equation of the above implies Substituting the above into the second equation of (3.6), we have Since ker − 1 2 I + (K −α,0 D ) * is generated by ψ α 1 and ψ α 2 , we may write φ as where |A| + |B| = 1. By integrating (3.8) on ∂D, and using the following identity Then, since we have from (2.7) that and ω * = (δc α So we get where, as in Theorem 2.1, Then, by solving the above eigenvalue problem, we obtain two (approximate) eigenpairs as follows: ω α ± = ω * ± λ δ ǫ|c| · |α| + O(δǫ 2 + δ 2 ), and (3.14) This implies that the Bloch eignfunction u α + (resp., u α − ) associated to the upper part ω α + (resp., ω α − ) of the Dirac cone can be represented as Then, by Corollary 3.2 and the fact that ω α

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. We see from the above lemma that the Dirac cone is located at the point α * /s. We denote the Dirac frequency as We have the following result for the Bloch eigenfunctions u α/s j,s , j = 1, 2 near the Dirac points α * /s.
Proof. The conclusion follows by applying (3.15) to u We see that the functions S 1 , S 2 describe the microscopic behaviour of the Bloch eigenfunction u α * /s+α ±,s while A ± e iα·x , B ± e iα·x describe the macroscopic behaviour. Now we derive a homogenized equation near the Dirac frequency ω * s . Recall that the Dirac frequency of the unscaled honeycomb crystal satisfies ω * ≈ C √ δ. As in [9], to make the order of ω * s fixed when s tends to zero, we assume that Assumption 4.1. δ = µs 2 , for some fixed µ > 0.
Then we have So, in what follows, we omit the subscript s in ω * s , namely, ω * := ω * s . Suppose the frequency ω is close to ω * , i.e., We need to find the Bloch eigenfunctions orα such that We have from (3.13) and Lemmas 4.1 and 4.2 that the correspondingα satisfies So it is immediate to see that the macroscopic field [ũ 1 ,ũ 2 ] T := [A ± e iα·x , B ± e iα·x ] T satisfies the Dirac equation as follows: Here, the superscript T denotes the transpose and ∂ i is the partial derivative with respect to the ith variable. Note that the each componentũ j , j = 1, 2 of the macroscopic field satisfies the Helmholtz equation The following is the main result on the homogenization theory for the honeycomb bubbly crystals. Theorem 4.3. For frequencies ω close to the Dirac frequency ω * , namely, ω−ω * = β √ δ, the following asymptotics of the Bloch eigenfunction u α * /s+α s holds: where the macroscopic field [ũ 1 ,ũ 2 ] T := [Ae iα·x , Be iα·x ] T satisfies the two-dimensional Dirac equation 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.

Numerical illustrations
In this section, we illustrate the main result, 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 also serves to illustrate and numerically verify the conclusions from [9]. The eigenmodes are computed by discretising the operator A α,ω δ from equation (2.11) 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 ρ = κ = 1000, ρ b = κ b = 1, which gives δ = 10 −3 , v = 1, and v b = 1.

Honeycomb lattice
For simplicity, we choose the scaling s = 1. Recall from (2.9) that the eigenmodes can be expressed as where A α,ω δ φ ψ = 0. Then (φ, ψ) can be numerically computed as an eigenvector of the discretised operator A α,ω δ corresponding to 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 ω with ω − ω * = β √ δ. First, we consider the small-scale behaviour of the eigenmodes. The small scale corresponds to e iα·x ≈ 1, so Theorem 4.3 shows that the eigenfunctions are given by Equation (3.14) shows that A = A ± = ± 1 √ 2 e i(θc+θ) and B = B ± = 1 √ 2 , where θ c , θ are the arguments of c andα respectively, and the sign coincides with the sign of β. To pick a unique eigenmode, we choose θ = 0, i.e., quasi-periodicityα = α 1 0 . Figure 3 shows the function which is the first Dirac eigenmode in the limit β → 0−. It can be seen that the eigenmode is highly oscillating and oscillate between −1 and 1 within one hexagon of bubbles. Moreover, this eigenmode has no large-scale oscillation. The second eigenmode, corresponding to β → 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 α * +α − for β = 8·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 iα·x S 1 (x) and Be iα·x S 2 (x). These fields are phase-shifted, due to the factor e i(θc+θ) in A.
To demonstrate the near-zero effective refractive index property of the bubbly honeycomb crystal, we consider the large-scale oscillation frequency close to the Dirac points. From (4.2), we know that for frequencies ω = ω * + β √ δ close to the Dirac frequency, each componentsũ j , j = 1, 2, of the macroscopic field oscillates at a spatial frequency Denote by ǫ := β √ δ. To verify this relation, we compute the large-scale spatial frequency of the eigenmodes for ǫ in the range ǫ ∈ [−0.01, 0.01], shown in Figure 5. It can be seen that the relation is linear for ǫ close to 0. This verifies (4.2), 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 Γ, 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.

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 ω * of the first band is achieved at the symmetry point α * = M = (π, π) in the Brillouin zone.
Let now S α,k D denote the single-layer potential for the square lattice, defined analogously as in Section 2.2 but with Λ and Λ * being the square lattice and reciprocal square lattice, respectively. Define the function    We now consider the rescaled crystal with unit cell sY, s > 0. Again, we assume that the order of ω * is fixed, i.e. δ = µs 2 , for some µ > 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 ω close to the critical frequency ω * , namely, (ω * ) 2 − ω 2 = O(s 2 ), the following asymptotics of the Bloch eigenfunction u α * /s+α s holds: where the macroscopic fieldũ := e iα·x satisfies the Helmholtz equatioñ which can be considered as a homogenized equation for the square bubbly crystal, while the microscopic field S vary on the scale of s.
The isotropic form of the macroscopic equation (5.2) follows since D is a circle, and an expression forλ is given in [9].
We now compute the eigenmodes of the square crystal close to the critical frequency, namely, ω = ω * − ǫ. The small-scale behaviour of the eigenmodes, i.e. the function S x s , 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α = α 1 0 and ω * − ω = 6 · 10 −5 . Similarly as in the case of a honeycomb, the eigenfunction varies at the large scale with a low-frequency macroscopic fieldũ.
To illustrate the macroscopic equation (5.2), the spatial frequency f of the macroscopic field was computed for ǫ in the range ǫ ∈ [0, 0.01]. Observe that due to the bandgap above ω * , no eigenmodes exists for ǫ < 0. Figure 9 shows the spatial frequency   f for different ǫ. This figure shows that f scales like √ ǫ for small ǫ > 0. This is consistent with (5.2), from which we expect a spatial frequency 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 behaviour close to the critical frequency, due to the fact that the square lattice cannot be mapped to a zero-index effective material.

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. Moreover, we have shown that the envelopes satisfy a Dirac equation. A comparison with a square lattice structure shows the great potential of using honeycomb crystals of sub-wavelength resonators as near-zero material. 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 Γ and therefore, may exhibit finite acoustic impedance and a high transmittance [14]. We also plan to mathematically analyse wave propagation phenomena in near-zero materials.