High-frequency homogenization of zero-frequency stop band photonic and phononic crystals

Over the last 25 years, many significant advances have created a deep understanding of the optical properties of photonic crystals (PCs) (John 1987, Yablonovitch 1987); such periodic structures prohibit the propagation of light, or allow it only in certain directions at certain frequencies, or localize light in specified areas. This sort of metamaterial (using the consensual terminology for artificial materials engineered to have desired properties that may not be found in nature, such as negative refraction, see e.g. Smith et al 2004, Ramakrishna 2005) enables a marked enhancement of control over light propagation; this arises from, for instance, the periodic patterning of small metallic inclusions embedded within a dielectric matrix (Nicorovici et al 1995, Pendry et al 1996). PCs are periodic devices whose spectrum is characterized by photonic band gaps and pass bands, just as electronic band gaps exist in semiconductors. In PCs, light propagation is disallowed for certain frequencies in certain (partial band gap, Yablonovitch and Gmitter 1989) or even all (complete band gap, Ho et al 1990) directions. This effect is well known, with many textbooks offering comprehensive introductions to the physics and models of PCs (Zolla et al 2005, Joannopoulos et al 2008, Markos and Soukoulis 2008, Ramakrishna and Grzegorczyk 2008), and forms the basis of many devices, including Bragg mirrors, dielectric Fabry–Perot filters and distributed feedback lasers. All of these contain low-loss dielectrics periodic in one dimension, so are one-dimensional PCs. Such mirrors are tremendously useful, but their reflecting properties critically depend upon the frequency of the incident wave with regard to its incidence (Markos and Soukoulis 2008). For broad frequency ranges, one wishes to reflect light of any polarization at any angle (which requires a complete photonic band gap) and for Dirichlet media (i.e. those composed with microstructure where the field is zero on the microparticles) such a gap occurs at zero frequency. A similar situation occurs in periodically perforated plates with clamped holes: so-called zero-frequency platonic stop bands disallow propagation of flexural waves at arbitrarily low frequencies (Antonakakis and Craster 2012). Using approximate models of thin plates to predict the propagation of surface elastic waves in structured soils (Brule et al 2012), it is therefore possible to create seismic shields for moderately low (Brule et al 2013) and even ultra-low-frequency waves (Antonakakis et al 2013b). Extending these ideas to higher dimensions allows for a much greater range of optical effects: all-angle negative refraction (Zengerle 1987, Luo et al 2002), ultra-refraction (Farhat et al 2010, Craster et al 2011) and cloaking at Dirac-like cones (Chan et al 2012). We will demonstrate here that an effective, frequency-dependent medium can be created that accurately reproduces these effects.

Considerable recent activity in this area is, in part, fuelled by advances in numerical techniques, based on Fourier expansions in the vector electromagnetic Maxwell equations (Joannopoulos et al 2008), finite elements (Zolla et al 2005) and multipole and lattice sums (Movchan et al 2002, Borwein et al 2013) for cylinders, to name but a few, that allows one to visualize the various effects and design PCs. The multipole method takes its root in a seminal paper of Lord Rayleigh (1892) which is the foundation stone upon which the current edifice of homogenization theories is built. More precisely, in that paper, John William Strutt, the third Lord Rayleigh, solved Laplace's equation in two dimensions for rectangular arrays of cylinders, and in three dimensions for cubic lattices of spheres. However, a limitation of Rayleigh's algorithm is that it does not hold when the volume fraction of inclusions increases. Multipole methods, in conjunction with lattice sums, overcome such obstacles and lead to the Rayleigh system which is an infinite linear algebraic system; this formulation, in terms of an eigenvalue problem, facilitates the construction of dispersion curves and the study of both photonic and phononic band-gap structures. In the limit of small inclusion radii, when propagating modes are very close to plane waves, one can truncate the Rayleigh system ignoring the effect of higher multipoles, to produce a series of approximations each successively more accurate to higher values of filling fraction. At the dipole order, one is able to fit the acoustic band, which has linear behaviour in the neighbourhood of zero frequency, except of course in the singular case whereby Dirichlet data are enforced on the boundary inclusions (Poulton et al 2001) and there is a zero-frequency stop band. This is a substantial limitation of this truncation. We will derive here the exact solution for zero-radius Dirichlet holes, avoiding multipole expansions, and show dispersion curves whose features lead to interesting topical effects, such as Dirac-like cone cloaking and degenerate points (Chan et al 2012), which emerge very naturally in this exact solution.

Although the numerical approaches discussed briefly are efficient, they still require substantial computational effort and can obscure physical understanding and interpretation. Here we substantially reduce the numerical complexity of the wave problem using asymptotic analysis; this has been developed over the last 35 years by applied mathematicians primarily for solving PDEs, with rapidly oscillating periodic coefficients, in the context of thermostatics, continuum mechanics or electrostatics (Bensoussan et al 1978, Bakhvalov and Panasenko 1989, Jikov et al 1994, Milton 2002). The available literature on such effective medium theories is vast, but it seems that only a very few groups have addressed such problems as the homogenization of media, with moderate contrast in the material properties, in three dimensions (Guenneau and Zolla 2000, Wellander and Kristensson 2003) and high-contrast two-dimensional (Zhikov 2000, Bouchitté and Felbacq 2004, Cherednichenko et al 2006) photonic crystals, that have important potential applications in photonics (e.g. microstructured fibres, Cregan et al 1999, Zolla et al 2005). Besides this, the aforementioned literature does not address the challenging problem of homogenization for moderate contrast photonic crystals near stop band frequencies. Classical homogenization is constrained to low frequencies and long waves in moderate contrast PCs (Silveirinha and Fernandes 2005) and so-called high-contrast homogenization only captures the essence of stop bands in PCs when the permittivity inside the inclusions is much higher than that of the surrounding matrix (the contrast being typically of the order of ⁻², where is the array pitch which in turn is much smaller than the wavelength). This latter area of homogenization theory is fuelled by interest in artificial magnetism, initiated by the work of O'Brien and Pendry (2002). However, moderate contrast one-dimensional photonic crystals have been recently shown to display not only artificial magnetism, but also chirality (Liu et al 2013), which is an onset to negative refraction (Pendry 2004).

For all these reasons, there is a strong demand for high-frequency homogenization (HFH) theories of PCs in order to grasp, and fully exploit, the rich behaviour of photonic band gap structures (Notomi 2000, Paul et al 2011). This has created a suite of extended homogenization theories for periodic media called Bloch homogenization (Conca et al 1995, Allaire and Piatnitski 2005, Birman and Suslina 2006, Hoefer and Weinstein 2011). There is also a flourishing literature on developing homogenized elastic media, with frequency-dependent effective parameters, also based upon periodic media (Nemat-Nasser et al 2011). There is considerable interest in creating effective continuum models of microstructured media that break free from the conventional low-frequency homogenization limitations.

In this paper, we apply the HFH theory proposed by one of us 3 years ago (Craster et al 2010) to the physics of zero-frequency band gap structures, not only in the context of photonics but also phononics for anti-plane shear waves in periodic arrays of inclusions, and platonics with flexural waves in pinned plates. The latter analysis is made possible thanks to the extension of HFH to plate theory, developed by two of us 2 years ago (Antonakakis and Craster 2012). Our aim is to show the universal features of stop band structures thanks to HFH, and to further exemplify their potential use in control of light and mechanical waves, with novel applications ranging from cloaking (Milton and Nicorovici 2006, Guenneau et al 2007) to anti-earthquake seismic shields (Antonakakis et al 2013a, Brule et al 2013).

We show in figure 1 what HFH does in practice. It focuses its attention on the physics within a supercell of sidelength L, which is the long scale, and further captures the fine features of field oscillations inside an elementary cell of sidelength 2l much smaller than L (Craster et al 2010). The wavelength need not be large compared to l in order to perform the asymptotic analysis, as the small parameter  = l/L only requires the supercell to be much larger than its constituent elementary cells; this contrasts with classical homogenization that assumes the cells are much smaller than the wavelength. In this way, one replaces a periodic structure by a homogenized (dispersive) one, on the long scale at any frequency, and classical homogenization is just a particular case of HFH (Antonakakis et al 2013a). An illustration of the HFH theory is in figure 2; the left panels are from full finite element simulations of an array of cylindrical holes (infinite conducting boundary conditions of the holes which are for transverse magnetic (TM) waves in electromagnetism, or clamped shear horizontal waves in elasticity) excited at the array centre by a line source and the right panels replace the array with its effective HFH medium which reproduces the large scale behaviour. For the frequencies chosen, two very distinct cross-shaped propagations, a Saint Andrews' cross (or saltire) (figures 2(a) and (b)) and a Saint George's cross (figures 2(c) and (d)) are achieved; these distinctive shapes are created by excitations at precise frequencies that are identified from the dispersion curves (Craster et al 2012), this strong frequency-dependent anisotropy is a recurrent feature of PCs. One sees that the HFH neatly reproduces in figures 2(b) and (d) the fine features of the full finite element simulations. Interpretation is provided by using the Brillouin zone of figure 3 and dispersion curves shown in figure 4. The square array of holes behaves as two dramatically contrasting effective media at frequency Ω = 1.966 (see figures 2(a) and (b)) which stands on the edge of the zero-frequency stop band that coincides with the lower edge of the second stop band; and at frequency Ω = 2.75 which stands on the upper edge of the second stop band, where the dispersion curve is nearly flat is the MX direction. Both of these frequencies are well beyond the range of applicability of classical homogenization.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our aim is to generate HFH for arrays of holes with Dirichlet conditions which represent TM waves in electromagnetism. This will enable us to compare, using full numerical simulations of HFH and finite elements, the behaviour of the numerous interesting optical effects which emerge through the coefficients of the effective medium. This is intertwined with knowledge and understanding of the wave structure for a perfect medium where the dispersion relations connecting the phase shift across an elementary cell with the frequency are important. In section 2 we set up the mathematical framework leading to the effective medium equation (16), this is in a dimensionless setting so for clarity in section 3 the main results are summarized in a dimensional setting. There are degenerate cases that are of substantial interest and connect with Dirac-like dispersion (Chan et al 2012) and lead to a different effective equation also outlined in section 3.

A first step in verifying the HFH theory is by asymptotically reproducing the dispersion curves for perfect infinite arrays. We use both circular and square holes for illustration, as shown in section 4. Near the edges of the Brillouin zone there are standing wave frequencies and the local asymptotic behaviour is given from the effective medium equation. Importantly, one can then predict a priori, from the signs of the HFH coefficients with the effective equation, how the bulk medium will behave. As a further illustration, holes that degenerate to a point are treated (section 4.3) as there is then an exact and simple solution for the dispersion relation; furthermore in this case degeneracies occur. Finally, in section 5 we compare HFH with full numerical simulations showing the full range of features available and how HFH represents them. Some concluding remarks are drawn together in section 6.

A two-dimensional structure composed of a doubly periodic, i.e. periodic in both x and y directions, square array of cells with, not necessarily circular, identical holes inside them is considered (see figure 3). We are primarily concerned with electromagnetism where Maxwell's equations separate into transverse electric (TE) and TM polarizations and the natural boundary conditions on the holes are then Neumann and Dirichlet, respectively; the asymptotics for the former are considered in Antonakakis et al (2013a) and applications in metamaterials are illustrated with split ring resonators (Pendry et al 1999). In the current paper our emphasis is rather different, and is upon TM polarization for which the Dirichlet conditions induce different effects, such as zero-frequency stop bands, and require a modified theory.

Assuming an $\exp (-\mathrm {i}\omega t)$ time dependence that is considered understood, and henceforth suppressed, the governing equation is

$$\begin{equation} \nabla_{{\bf x}} [\hat{a}({\boldsymbol{\xi}})\nabla_{{\bf x}}u({\bf x})]+\omega^2\hat{\rho}({\boldsymbol{\xi}})u({\bf x})=0 \end{equation} \tag{ 1 }$$

on −∞ < x₁, x₂ < ∞, where ω is the angular frequency. In the periodic setting, each cell is identical and the material is characterized by two periodic functions on the short scale, in ξ ≡ (x₁/l,x₂/l), namely $\hat {a}({\boldsymbol{\xi}})$ and $\hat {\rho }({\boldsymbol{\xi}})$. In the context of photonics, these are the inverse of permeability (μ⁻¹) and permittivity (ε), respectively, which are related to the wavespeed of light in vacuum c via the refractive index n as follows: εμ = n²c⁻². Note that the permittivity (ε) is not to be confused with the small parameter . The unknown u physically is the longitudinal component E_z of the out-of-plane electric field E = (0,0,E_z), bearing in mind that the transverse magnetic field H = (H_x,H_y,0) is retrieved from Maxwell's equation H = iω⁻¹μ⁻¹∇ ×  E. In the context of phononics, $\hat {a}({\boldsymbol{\xi}})$ and $\hat {\rho }({\boldsymbol{\xi}})$ would stand, respectively, for the shear modulus and the density of an isotropic elastic medium, and the unknown u would correspond to the out-of-plane displacement (amplitude of SH shear waves). One can also easily draw analogies with acoustic pressure waves in a fluid, or linear surface water waves, which explains the activity related to finding correspondences between models of electromagnetic (Ramakrishna 2005) and acoustic (Craster and Guenneau 2013) metamaterials. The lengths of the direct lattice base vectors are taken equal to 2l, and define a short length scale. The overall dimension of the structure is of a much longer length scale, L; the ratio of scales,  ≡ l/L ≪ 1 then provides a small parameter for use later in the asymptotic scheme.

A non-dimensionalization of the physical functions by setting $\hat {a}\equiv \hat {a}_{0}a({\boldsymbol{\xi}})$ and $\hat {\rho }\equiv \hat {\rho }_{0}\rho ({\boldsymbol{\xi}})$ is convenient. Equation (1) then becomes

$$\begin{equation} l^2\nabla_{{\bf x}} [a({\boldsymbol{\xi}})\nabla_{{\bf x}}u({\bf x})]+\Omega^2\rho({\boldsymbol{\xi}})u({\bf x})=0 \quad {\rm with} \; \Omega = \frac{\omega l}{\hat{c}_{0}} \end{equation} \tag{ 2 }$$

as the dimensionless frequency and $\hat {c}_0=\sqrt {\hat {a}_0/\hat {\rho }_0}$. The two scales, l,L, then motivate two new coordinates, namely X =  x/L and ξ =  x/l that are treated as independent, and placing this change of coordinates into (2) we obtain

$$\begin{eqnarray} &&\fl \nabla_{{\boldsymbol{\xi}}}\cdot[a({\boldsymbol{\xi}})\nabla_{{\boldsymbol{\xi}}}u({\bf X},{\boldsymbol{\xi}})]+\Omega^2\rho({\boldsymbol{\xi}}) u ({\bf X},{\boldsymbol{\xi}}) +\epsilon[2a({\boldsymbol{\xi}})\nabla_{\boldsymbol{\xi}} +\nabla_{{\boldsymbol{\xi}}}a({\boldsymbol{\xi}})]\cdot\nabla_{{\bf X}}u({\bf X},{\boldsymbol{\xi}})\nonumber\\ &&+\epsilon^2 a({\boldsymbol{\xi}})\nabla^2_{\bf X}u({\bf X},{\boldsymbol{\xi}})=0. \end{eqnarray} \tag{ 3 }$$

For wave propagation through a perfect lattice there exist standing waves at specific eigenfrequencies, Ω₀, and their eigenmodes satisfy periodic or anti-periodic (out-of-phase) boundary conditions on the edges of the cell. We obtain asymptotic solutions based upon these standing waves.

A key point, from this separation of scales, is that the boundary conditions on the short scale are known. For instance, periodic conditions in ξ on the edges of the cell are

$$\begin{equation} u|_{\xi_{i}=1}=u|_{\xi_{i}=-1} \quad {\rm and} \quad u_{,\xi_i}|_{\xi_{i}=1}=u_{,\xi_i}|_{\xi_{i}=-1} , \end{equation} \tag{ 4 }$$

where u_,ξi denotes partial differentiation with respect to rescaled coordinate ξ_i. An almost identical analysis holds near the anti-periodic standing wave frequencies. We illustrate the periodic case for definiteness. On the long scale we set no boundary conditions, and indeed the methodology we present can be used even when the problem is not perfectly periodic (Craster et al 2011).

The following ansatz is taken:

$$\begin{equation} \fl u({\bf X},{\boldsymbol{\xi}})=u_0({\bf X},{\boldsymbol{\xi}})+\epsilon u_1({\bf X},{\boldsymbol{\xi}})+\epsilon^2 u_2({\bf X},{\boldsymbol{\xi}})+\cdots , \quad \Omega^2=\Omega_0^2+\epsilon \Omega_1^2+\epsilon^2 \Omega_2^2+\cdots \end{equation} \tag{ 5 }$$

Since u( X,ξ) is periodic in ξ so are the components u_i( X,ξ)'s. This leads to a hierarchy of equations in increasing powers of : the leading order followed by the first order, second order and so on. The first three equations in ascending order yield

$$\begin{equation} (au_{0,\xi_i})_{,\xi_i} + \Omega_0^2\rho u_0=0, \end{equation} \tag{ 6 }$$

$$\begin{equation} (au_{1,\xi_i})_{,\xi_i} + \Omega_0^2\rho u_1= -(2au_{0,\xi_i}+a_{,\xi_i}u_0)_{,X_i} -\Omega_1^2 \rho u_0, \end{equation} \tag{ 7 }$$

$$\begin{eqnarray} &&(au_{2,\xi_i})_{,\xi_i} + \Omega_0^2\rho u_2 =-au_{0,X_iX_i} -(2au_{1,\xi_i}+ a_{,\xi_i}u_1)_{,X_i}-\Omega_1^2\rho u_1 -\Omega_2^2\rho u_0. \end{eqnarray} \tag{ 8 }$$

We start by expressing a solution for the leading order equation. For a specific eigenvalue Ω₀ there is a corresponding eigenmode U₀(ξ;Ω₀), periodic in ξ, and the leading order solution is

$$\begin{equation} u_0({\bf X},{\boldsymbol{\xi}})=f_0({\bf X}) U_0({\boldsymbol{\xi}};\Omega_0). \end{equation} \tag{ 9 }$$

The function f₀( X) is determined later on, indeed the whole point is to deduce a long-scale continuum partial differential equation for f₀ entirely upon the long-scale X. We begin with a treatment of isolated eigenvalues, however repeated eigenvalues also arise, and require subtle modifications, and will be dealt with in sections 2.2 and 2.3. These are relevant to Dirac-like cones (Chan et al 2012).

Dirichlet conditions on the inside boundary of the cell, ∂S₂, where S denotes the surface of the cell in the ξ coordinates, yield

$$\begin{equation} u({\bf X},{\boldsymbol{\xi}})|_{\partial S_2}=0 \iff u_i({\bf X},{\boldsymbol{\xi}})|_{\partial S_2}=0, \quad i\in\ \mathbb{N} \end{equation} \tag{ 10 }$$

and are set in the short-scale ξ so, for i = 0, U₀(ξ;Ω₀)|_∂S2 = 0.

To deduce the long-scale PDE for f₀( X) we follow Craster et al (2010), making changes where necessary due to the change in boundary conditions. Equation (7) is multiplied by U₀ and integrated over the cell's surface to obtain

$$\begin{eqnarray} &&\fl \int\!\!\!\!\int_S \left( U_0(a u_{1,\xi_i})_{,\xi_i} +\Omega_0^2\rho U_0 u_1\right) \mathrm{d}S =-f_{0,X_i} \int\!\!\!\!\int_S(a U_0^2)_{,\xi_i}\,\mathrm{d}S -f_0\Omega_1^2 \int\!\!\!\!\int_S\rho U_0^2\, \mathrm{d}S . \end{eqnarray} \tag{ 11 }$$

The first integral of the right-hand side is converted to a path integral along ∂S using a corollary of Green's theorem. The periodic conditions on ∂S₁ and the homogeneous Dirichlet conditions of U₀ on ∂S₂ make it vanish. We continue by subtracting the integral over the cell of the product of equation (6) with u₁/f₀ to obtain

$$\begin{equation} \int\!\!\!\!\int_S (U_0(a u_{1,\xi_i})_{,\xi_i} - u_1(a U_{0,\xi_i})_{,\xi_i})\, \mathrm{d}S = -f_0\Omega_1^2 \int\!\!\!\!\int_S \rho({\boldsymbol{\xi}})U_0^2 \,\mathrm{d}S. \end{equation} \tag{ 12 }$$

Using Green's theorem equation (12) becomes

$$\begin{equation} \int_{\partial S} a({\boldsymbol{\xi}})\left(U_0\frac{\partial u_1}{\partial n}-u_1\frac{\partial U_0}{\partial n}\right)\, \mathrm{d}s=-f_0\Omega_1^2 \int\!\!\!\!\int_S \rho({\boldsymbol{\xi}})U_0^2\, \mathrm{d}S, \end{equation} \tag{ 13 }$$

where ∂S = ∂S₁∪∂S₂ and ∂/∂n: =  n∇ with n the outer normal to ∂S. Due to periodic (or anti-periodic) boundary conditions of u and its first derivatives with respect to ξ_i on ∂S₁, and due to the homogeneous Dirichlet boundary conditions of the u_i's on ∂S₂, the left-hand side of equation (13) vanishes. It follows that Ω₁ = 0, which is an important deduction as it implies that the asymptotic behaviour of dispersion curves near isolated eigenfrequencies is at least quadratic. The right hand side of equation (7) reduces to −(2au_0,ξi + a_,ξiu₀)_,Xi which implies the presence of a f_0,Xi term in the solution for u₁( X,ξ) that reads

$$\begin{equation} u_1({\bf X},{\boldsymbol{\xi}})=f_1({\bf X}) U_0({\boldsymbol{\xi}};\Omega_0) +\nabla_{\bf X} f_0({\bf X})\cdot {\bf U_1}({\boldsymbol{\xi}}). \end{equation} \tag{ 14 }$$

The first term is simply a homogeneous solution that is irrelevant, and the vector function U₁ is a particular solution. Inserting equation (14) into equation (7) the solution of U_1j reads

$$\begin{equation} (aU_{1j,\xi_i})_{,\xi_i} + \Omega_0^2\rho U_{1j}= -2aU_{0,\xi_j}-a_{,\xi_j}U_0\quad{\rm for}\; j=1,2. \end{equation} \tag{ 15 }$$

Quite often exact solutions for u₁ or U₁ are not possible to obtain so one must use numerical solutions and this is described in section 5. We now move to the second order equation to compute f₀( X) and Ω²₂. Similar to the first order equation, we multiply equation (8) by U₀ then subtract the product of equation (7) with u₂/f₀, finally we integrate over the cell's surface to obtain the following effective equation for f₀:

$$\begin{equation} T_{ij} f_{0,X_iX_j}+\Omega_2^2 f_0=0,\quad {\rm with}\; T_{ij}=\frac{t_{ij}}{\int\!\!\int_S \rho U_0^2\, \mathrm{d}S}\quad {\rm for} \quad i,j=1,2, \end{equation} \tag{ 16 }$$

which is the main homogenization result of this paper. In particular, the coefficients T_ij encode the anisotropy of the effective medium at a specific frequency and the t_ij's are integrals over the small-scale cell

$$\begin{equation} t_{ii}=\int\!\!\!\!\int_SaU_0^2\,\mathrm{d}S+2\int\!\!\!\!\int_SaU_{1_i,\xi_i}U_0\,\mathrm{d}S+\int\!\!\!\!\int_S a_{,\xi_i}U_{1_i}U_0\,\mathrm{d}S, \end{equation} \tag{ 17 }$$

$$\begin{equation} t_{ij}=2\int\!\!\!\!\int_SaU_{1_j,\xi_i}U_0\,\mathrm{d}S+\int\!\!\!\!\int_S a_{,\xi_i}U_{1_j}U_0\,\mathrm{d}S \quad {\rm for}\quad i\neq j, \end{equation} \tag{ 18 }$$

where U_1i is the ith component of vector function U₁. Note that there is no summation over repeated indices for t_ii. The PDE for f₀, equation (16), is crucial as the local microstructure is completely encapsulated within the tensor T_ij; these are, for a specific structure at an Ω₀, just numerical values as illustrated in table 2. Notably the tensor can have negative values, or components, and it allows one to interpret and, even more importantly, predict changes in behaviour when the frequency changes, or when specific effects occur. The structure of the tensor depends upon the boundary conditions of the holes: the results here for instance, are markedly different from those of the Neumann case (Antonakakis et al 2013a). Another key point is that numerically the short scale is no longer present and the PDE (16) is simple and quick to solve numerically, or even by hand.

One primary aim of the present paper is to deduce formulae for the local behaviour of dispersion curves near standing wave frequencies as these shed light upon the physical effects observed. For this, one returns to the perfect lattice and then Floquet–Bloch boundary conditions on the elementary cell lead immediately to $f_0({\bf{ X}})=\exp (\mathrm {i}\kappa _j X_j/\epsilon )$. In this notation κ_j = K_j − d_j and d_j = 0, π/2, −π/2 depending on the location in the Brillouin zone, see figure 3(b). Equation (16) gives

$$\begin{equation} \Omega\sim\Omega_0+ \frac{T_{ij}}{2\Omega_0}\kappa_i \kappa_j \end{equation} \tag{ 19 }$$

and these locally quadratic dispersion curves are completely described by Ω₀ and the tensor T_ij.

2.1. Classical singularly perturbed zero-frequency limit

2.2. Repeated eigenvalues: linear asymptotics

Repeated eigenvalues are commonplace in specific examples and, as we shall see, are particularly relevant to Dirac-like cones (Chan et al 2012) that have practical significance. If we assume repeated eigenvalues of multiplicity p the general solution for the leading order problem becomes

$$\begin{equation} u_0=f_0^{(l)}({\bf X})U_0^{(l)}({\boldsymbol{\xi}};\Omega_0), \end{equation} \tag{ 20 }$$

where we sum over the repeated superscripts (l). Proceeding as before, we multiply equation (7) by U^(m)₀, subtract u₁((aU^(m)_0,ξi)_,gξi + Ω²₀ρU^(m)₀), then integrate over the cell to obtain

$$\begin{equation} f_{0,X_i}^{(l)}\int\!\!\!\!\int_SU_0^{(m)}\left(2aU_{0,\xi_i}^{(l)}+a_{\xi_i}U_0^{(l)}\right)\,\mathrm{d}S+\Omega_1^2 f_0^{(l)}\int\!\!\!\!\int_S\rho U_0^{(l)}U_0^{(m)}\,\mathrm{d}S=0. \end{equation} \tag{ 21 }$$

There is now a system of coupled PDEs for the f^(l)₀ and, provided Ω₁ ≠ 0, the leading order behaviour of the dispersion curves near the Ω₀ is now linear (these then form Dirac-like cones). These coupled partial differential equations on the long scale now replace (16) near these frequencies.

For the perfect lattice and Bloch wave problem, we set $f_0^{(l)}={\skew5\hat f}_0^{(l)}\exp (\mathrm {i}\kappa _jX_j/\epsilon )$ and obtain the following equations:

$$\begin{equation} \left(\mathrm{i}\frac{\kappa_j}{\epsilon}{\bf A}_{jml}+\Omega_1^2{\bf B}_{ml}\right){\skew5\hat{f}}_{\!0}^{(l)}=0, \quad {\rm for} \quad m=1,2,\ldots,p, \end{equation} \tag{ 22 }$$

where

$$\begin{equation} {\bf A}_{jml}=\int\!\!\!\!\int_SU_0^{(m)}\left(2aU_{0,\xi_j}^{(l)}+a_{,\xi_j}U_0^{(l)}\right)\,\mathrm{d}S \end{equation} \tag{ 23 }$$

and

$$\begin{equation} {\bf B}_{ml}=\int\!\!\!\!\int_S\rho U_0^{(l)}U_0^{(m)}\,\mathrm{d}S. \end{equation} \tag{ 24 }$$

The system of equations (22) is written simply as

$$\begin{equation} {\bf C}{\hat {\bf F}}_0={\bf 0}, \end{equation} \tag{ 25 }$$

with C_ll = Ω²₁ B_ll and C_ml = iκ_j A_jml/ for l ≠ m. One must then solve for $\Omega _1^2=\pm \sqrt {\alpha _{ij}\kappa _i\kappa _j}/\epsilon$ when the determinant of C vanishes and the asymptotic relation is

$$\begin{equation} \Omega\sim\Omega_0\pm\frac{1}{2\Omega_0}\sqrt{\alpha_{ij}\kappa_i\kappa_j}. \end{equation} \tag{ 26 }$$

If Ω₁ is zero, one must go to the next order and a slightly different analysis ensues.

2.3. Repeated eigenvalues: quadratic asymptotics

Assuming that Ω₁ is zero, u₁ = f^(l)_0,XkU^(l)_1k (we again sum over all repeated (l) superscripts) and advance to second order using (8). Taking the difference between the product of equation (8) with U^(m)₀ and u₂((aU_0,ξi)_,ξi + Ω²₀ρU₀) and then integrating over the elementary cell gives

$$\begin{eqnarray} &&\fl f_{0,X_iX_i}^{(l)}\int\!\!\!\!\int_SaU_0^{(m)}U_0^{(l)}\,\mathrm{d}S+f_{0,X_kX_j}^{(l)}\int\!\!\!\!\int_SU_0^{(m)}\left(2aU_{1_k,\xi_j}^{(l)}+a_{,\xi_j}U_{1_k}^{(l)}\right)\,\mathrm{d}S\nonumber\\ &&+\Omega_2^2 f_0^{(l)}\int\!\!\!\!\int_S\rho U_0^{(m)}U_0^{(l)}\,\mathrm{d}S=0, \quad {\rm for} \quad m=1,2,\ldots,p \end{eqnarray} \tag{ 27 }$$

as a system of coupled PDEs. The above equation is presented more neatly as

$$\begin{equation} f_{0,X_iX_i}^{(l)}{\bf A}_{ml}+f_{0,X_kX_j}^{(l)}{\bf D}_{kjml}+\Omega_2^2 f_0^{(l)}{\bf B}_{ml}=0, \quad {\rm for} \quad m=1,2,\ldots,p. \end{equation} \tag{ 28 }$$

For the Bloch wave setting, using $f_0^{(l)}({\bf{ X}})={\skew5\hat{f}}_{\!0}^{(l)}\exp (\mathrm {i}\kappa _jX_j/\epsilon )$ we obtain the following system:

$$\begin{equation} \left(-\frac{\kappa_i\kappa_i}{\epsilon^2}{\bf A}_{ml}-\frac{\kappa_k\kappa_j}{\epsilon^2}{\bf D}_{kjml}+\Omega_2^2{\bf B}_{ml}\right){\skew5\hat{f}}_{\!0}^{(l)}=0, \quad {\rm for}\quad m=1,2,\ldots,p \end{equation} \tag{ 29 }$$

and this determines the asymptotic dispersion curves.

The major goal of HFH is to represent the periodic medium of finite extent using an effective homogeneous medium and the main theoretical result of this paper is in achieving this. For clarity this is summarized in this section back in the dimensional setting. An important assumption when deriving equation (16) is that the medium is perfectly periodic and of infinite extent. We will use equation (16) not only to homogenize infinite media (spectral problem) but also finite size periodic media embedded in a bulk material of infinite extent (scattering problem), and show that it represents a very good approximation. Transforming equation (16) back to the original x_i = X_iL coordinates and using the solution of Ω₂, we obtain an effective medium equation for ${\skew5\hat{f}}_{\!0}$ that is

$$\begin{equation} T_{ij}{\skew5\hat{f}}_{0,x_ix_j}({\bf x})+\frac{\Omega^2-\Omega_0^2}{l^2}{\skew5\hat{f}}_{\!0}({\bf x})=0. \end{equation} \tag{ 30 }$$

The nature of equation (30) is not necessarily elliptic since T₁₁ and T₂₂ can take different values and/or different signs. In the illustrations herein the T_ij = 0 for i ≠ j. The hyperbolic behaviour of equation (30) yields asymptotic solutions that describe the endoscope effects where, if one of the T_ii coefficients is zero as often happens near the point M of the Brillouin zone, waves propagate only in one direction. Note that when the cell has adequate symmetries the dispersion relation near point G(0,π/2) is identical to that near point X(π/2,0), which explains the existence of two orthogonal directions of propagation instead of only one. One can rewrite equation (30) in the form of

$$\begin{equation} C_{ij}{\skew5\hat{f}}_{0,x_ix_j}({\bf x})+\Omega^2{\skew5\hat f}_0({\bf x})=0, \end{equation} \tag{ 31 }$$

where C_ij = T_ijl²Ω²/(Ω² − Ω²₀) and can be identified as an effective product of inverse permittivity and inverse permeability.

At Dirac-like cones the linear behaviour of the effective medium yields an equation slightly different from (30). Regarding the zero radius holes presented in following section 4.3, this equation consists of a system of three coupled PDEs that uncouple at Dirac-like cones to yield one identical PDE for all f⁽ⁱ⁾₀s that is of the form

$$\begin{equation} {\skew5\hat{f}}_{0,x_ix_i}({\bf x})+\beta\frac{(\Omega^2-\Omega_0^2)^2}{l^2}{\skew5\hat{f}}_{\!0}({\bf x})=0, \end{equation} \tag{ 32 }$$

where β is a coefficient equivalent to the T_ij but this time is the same for all combinations of i and j. The quadratic mode that emerges in the middle of the linear ones follows equation (30). By rewriting equation (32) an effective refractive index emerges as in Craster et al (2011) that reads

$$\begin{equation} {\skew5\hat{f}}_{0,x_ix_i}({\bf x})+\eta_{\rm eff}\Omega^2{\skew5\hat{f}}_{\!0}({\bf x})=0, \end{equation} \tag{ 33 }$$

where

$$\begin{equation} \eta_{\rm eff}=\beta(\Omega^2-\Omega_0^2)^2/(\Omega^2l^2). \end{equation} \tag{ 34 }$$

Note that η_eff can take vanishing, or negative, values for frequencies just above or below Ω₀.

The effective medium parameters are given for wavenumbers that perturb away from the edges of the Brillouin zone. These parameters will be asymptotically exact as the excitation frequency Ω tends to the standing wave frequency Ω₀. HFH has no upper-limit frequency for Ω₀, although one should be aware of the numerical issues that can arise when solving for the standing wave eigenmodes: space discretization should be refined as Ω₀ reaches relatively high values in order to capture the finer structure of U₀ and U₁.

We now generate dispersion curves for circular and square holes in a matrix with properties a = ρ = 1 and verify the HFH theory versus these numerically generated curves; this is a good test of the approach as the dispersion curves contain changes in curvature, coalescing branches and modes that cross. These dispersion curves can then be used to interpret the physical optics phenomena seen later. Most vividly the dispersion curves also show the zero frequency, and other, stop bands.

4.1. Circular inclusions in a square array

The circular inclusions are arranged in a square array and each elementary cell is a square of side 2. We consider large holes initially and then how the curves change as the radius decreases. Ultimately we consider infinitesimal radii and perhaps remarkably obtain explicit solutions. In almost all other cases one has to use fully numerical techniques, or semi-analytical methods such as multipoles (Movchan et al 2002, Zolla et al 2005).

4.1.1. Large circular holes

The full dispersion curves are computed numerically using COMSOL (finite element code) and then asymptotically with the HFH using the equations developed in section 2. Figure 4 illustrates the typical dispersion curves versus the asymptotics for wavenumbers, on the specified path of figure 3, for the example of a hole with radius 0.4. There is, as expected, a zero-frequency stop band and the lowest branch is isolated with a full stop band also lying above it. At some standing wave frequencies two modes share the same frequency, for instance the third and fourth modes at Γ; these modes are asymptotically quadratic in the local wavenumber. Table 2 shows the T₁₁,T₂₂ values for point Γ of the Brillouin zone. Note how T₁₁ = T₂₂ for all single eigenfrequencies, but that symmetry breaks for the double roots. Moreover the signs of T₁₁ and T₂₂ naturally inform one of the local curvature near Γ. Physically, this tells us the sign of group velocity of waves with small phase-shift across the unit cells, and thus whether or not they undergo backward propagation, which is one of the hallmarks of negative refraction (Ramakrishna 2005).

4.1.2. Folding technique for reconstructing the dispersion curves asymptotically

An apparent deficiency of the approach we present is that it requires known standing wave frequencies, and associated eigenstates, at the band edges and that the asymptotics may be poor at frequencies far from these standing waves; until now HFH, applied to the dispersion curves, has been limited to obtaining asymptotic solutions near the standing wave frequencies (Craster et al 2010). This uses an elementary cell containing a single hole and the irreducible Brillouin zone associated with it; however, by using larger macro-cells containing four or more holes, and foldings of their resultant Brillouin zone, one can extract the asymptotics at other positions in wavenumber space thereby considerably enhancing the asymptotic coverage (see figure 5). As an example, we take an elementary square cell of length 2 with circular holes of radius 0.4, section 4.1 and figure 4.

Zoom In Zoom Out Reset image size

To obtain asymptotic solutions at the wavenumber positions, I(π/4,0), J(π/2,π/4), M'(π/4,π/4) of figure 5(b), bisecting the three segments [ΓX], [MX] and [ΓM] of the Brillouin zone we can use equations (16)–(18) again (section 2.2) but now with different elementary cells.

Let the Brillouin zone for a square macro-cell composed of four holes be Γ'M'I. At Γ'(0,0) periodic (in-phase) conditions on the macro-cell boundaries yield all in-phase and out-of-phase modes of the elementary cell containing a single hole. This elementary cell's [ΓM] segment is folded on point M'. At M', where anti-periodic (out-of-phase) conditions on the macro-cell are applied, the dispersion curves contain repeated roots, which are linear, and their slope is equal to the asymptote of the dispersion curves of the elementary cell setting at point M'. Proceeding in a similar manner one can fold the Brillouin zone at will by considering macro-cells composed of simple unfoldings of the elementary cell. Asymptotics can be obtained for wavenumber positions of the Brillouin zone equal to all ratios of the elementary cell's Brillouin zone wavenumber points.

We now illustrate this by obtaining asymptotics at points M', I and J which, respectively, correspond to the following macro-cell settings, a square macro-cell composed of four elementary cells with anti-periodic boundary conditions, a rectangular macro-cell composed of two elementary cells in the ξ₁ direction with anti-periodic boundary conditions on ξ₁ but periodic on ξ₂ and finally a rectangular macro-cell composed of two elementary cells in the ξ₂ direction with anti-periodic conditions in ξ₁ and ξ₂. The asymptotics at points M', I, J as well as the usual points Γ, M and X are illustrated in figure 5.

4.2. Circular inclusions of small radius

As the hole radius decreases, the dispersion curves gradually shift towards those of the infinitesimal holes shown in figure 6: the zero-frequency stop band is generic, and a feature forced by the Dirichlet boundary condition; however, the lowest dispersion curve is no longer isolated and triple crossings created by repeated roots occur.

Zoom In Zoom Out Reset image size

For the larger hole radius (section 4.1) all the modes are locally quadratic close to the standing wave frequencies. As the wavenumber moves away from those specific Brillouin zone points the behaviour of the dispersion curves becomes locally linear. The formation of triple, and more, crossings as the radius tends to zero, are created by the linear behaviour of the dispersion curves far away from points Γ, M and X moving towards those points to replace the quadratic behaviour of all the multiple modes, except one, as seen in figure 6. Triple crossings are illustrated at the second and fourth standing wave frequencies at point Γ and the first standing wave frequencies at point M, where the third frequency at point M yields a heptuple crossing. Reassuringly, the asymptotic theory captures these multiple modes whether the crossings are double, triple or more, as is illustrated in figures 4 and 6 where the circular holes are, respectively, of radius 0.4 and 0 (constrained points). For the multiple crossing points with multiplicity higher than two, both linear and quadratic dispersion curves coexist. One can proceed, as in section 2.2, to obtain the different Ω₁ coefficients and one of them is zero; proceeding to higher order and applying equations (27)–(29) gives the quadratic behaviour of the middle mode emerging from the triple crossings. Again the T_ii contain critical information about the local behaviour; some typical values along Γ are given in table 3, for the quadratic curves. They are naturally identical for isolated modes but describe the inherent anisotropy for the other cases. One can then immediately see whether the material will behave with directional preferences at specific frequencies.

Dispersion diagrams of arrays of constrained points and 0.1 radius holes (not shown) are almost indistinguishable by naked eye, but apparent triple crossings in the latter are actually a double crossing together with a single mode that has almost the same standing wave frequency. The separation of the two nearly touching standing wave frequencies depends on the size of the inclusion and tends to zero as the inclusion shrinks to a point. Figure 7 illustrates a close-up of a triple point comparing arrays of constrained points with inclusions of radius 0.1. This is actually important for the Dirac-like cones and associated effects.

Zoom In Zoom Out Reset image size

4.3. Constrained points, Fourier series and Dirac-like cones

Similar to arrays of constrained points for plates (Mace 1996), a doubly periodic array of constrained points (limit of circles with radius tending to zero) is an attractive limit and has an immediate solution in terms of Fourier series as

$$\begin{equation} u(\xi_1,\xi_2)=\mathrm{e}^{\mathrm{i}{\boldsymbol{\kappa}}\cdot{\boldsymbol{\xi}}}\sum_{n_1=-\infty}^\infty\sum_{n_2=-\infty}^\infty \frac{\mathrm{e}^{-\mathrm{i}\pi {\bf N}\cdot{\boldsymbol{\xi}}}}{{{\mid}{\boldsymbol \kappa}-\pi {\bf N}{\mid}}^2-\Omega^2} \end{equation} \tag{ 35 }$$

with κ = (κ₁,κ₂) and N = (n₁,n₂) and Ω and κ are related via the dispersion relation

$$\begin{equation} D({\boldsymbol{\kappa}},\Omega)=\sum_{n_1=-\infty}^\infty\sum_{n_2=-\infty}^\infty \frac{1}{{{\mid}{\boldsymbol \kappa}-\pi {\bf N}{\mid}}^2-\Omega^2}=0. \end{equation} \tag{ 36 }$$

This is valid provided the double sum does not have a singularity, which would occur whenever

$$\begin{equation} \Omega=\sqrt{{{\mid}{\boldsymbol \kappa}-\pi{\bf N}{\mid}}^2}, \end{equation} \tag{ 37 }$$

where |κ − π N|² = (κ₁ − πn₁)² + (κ₂ − πn₂)². These relations from the singularities also play a role in the dispersion picture; they correspond to the Bloch states of a perfect medium without any constraints (a homogeneous isotropic medium). However, in some circumstances they also, by serendipity, exactly satisfy the constraint u = 0 at the centre of each cell and therefore, in those cases, are simultaneously dispersion curves; this is the origin of the degeneracy seen in the multiple crossings (occurring at Dirac-like points) of the dispersion curves of figure 6. It is also clear from the dispersion relation that Ω = 0 is not a solution at κ =  0 and hence there is a zero-frequency stop band and conventional long-wave homogenization is doomed to fail (cf section 2.1).

A set of dispersion curves is shown in figure 6. The Bloch states for the perfect medium free of defects lead to light lines folded at the edges of the Brillouin zone and those emerging from Γ are pertinent to the current discussion. These light lines intersect at the multiple crossings, which we call generalized Dirac-like points. These occur for a given set of frequencies at all three edges of the Brillouin zone, namely for points Γ, M and X. It is useful to find criteria for their multiplicity. At Γ the generalized Dirac-like points occur for Ω² = π²m such that m = n²₁ + n²₂ with m integer; the multiplicity depends on the ways in which n₁ and n₂ can form m. From elementary number theory m can be equal to the sum of two squares in more than one way, i.e. m = 50 = 5² + 5² = 7² + 1² or the sum of two squares can also be a perfect square, i.e. m = 25 = 5² = 4² + 3². Equation u_,xixi + π²(n²₁ + n²₂)u = 0 together with the appropriate boundary conditions describes the generalized Dirac-like points and an independent set of solutions is formed from different possibilities of sin(π(n_ix + n_jy)), sin(π(n_ix − n_jy)) and cos(π(n_ix + n_jy)) − cos(π(n_px − n_ky)), where i,j,p,k∈{1,2}, i ≠ j and p ≠ k. The multiplicity can be equal to 3, 7, 11 or more depending on m. For example, Ω² = π² with n_i = 1, n_j = 0 has multiplicity 3 where Ω² = 5π² has multiplicity 7 and Ω² = 50π² has multiplicity 11. This is true for solutions at Γ and M but at X one can obtain a singularity that accepts only one eigensolution. Due to the non-symmetry of the general sum at X that renders the denominator of the Fourier series singular, m = (1/2 − n₁)² + n²₂, it is possible to obtain solutions of multiplicity one as in the case of Ω₀ = π/2 where the only eigensolution respecting the boundary conditions is sin(πx/2).

Given the exact solution one can, in these special cases, generate the asymptotics by hand. For the asymptotics we construct the coefficients using U₀ from (35) and deduce U₁ = (U₁₁,U₁₂) as

$$\begin{equation} U_{1k}=2\mathrm{i}\,\mathrm{e}^{\mathrm{i}{\boldsymbol{\kappa}}\cdot{\boldsymbol{\xi}}} \sum_{n_1=-\infty}^\infty\sum_{n_2=-\infty}^\infty \frac{(\kappa_k-\pi n_k)\, \mathrm{e}^{-\mathrm{i}\pi {\bf N}\cdot{\boldsymbol{\xi}}}}{[{{\mid}{\boldsymbol \kappa}-\pi{\bf N}{\mid}}^2-\Omega^2]^2} \end{equation} \tag{ 38 }$$

and thus one can extract the T_ij in (16) by doing the integrals by hand, the off-diagonal terms are zero, and T₁₁ and T₂₂ are

$$\begin{equation} T_{11}=1-4\frac{\sum_{n_1=-\infty}^\infty\sum_{n_2=-\infty}^\infty \frac{(\kappa_1-\pi n_1)^2}{[{\mid{\boldsymbol \kappa}-\pi{\bf N}\mid}^2-\Omega_0^2]^3}}{\sum_{n_1=-\infty}^\infty\sum_{n_2=-\infty}^\infty\frac{1}{[{\mid{\boldsymbol \kappa}-\pi{\bf N}\mid}^2-\Omega_0^2]^2}}, \end{equation} \tag{ 39 }$$

and an identical equation for T₂₂ but with 1 and 2 interchanged, where (κ₁,κ₂) are (0,0), (π/2,0) and (π/2,π/2) at the edges of the Brillouin zone at the respective points of Γ, X and M. These asymptotics only apply at the isolated, non-repeating roots, and are shown in figure 6.

The linear cases at the triple crossings are easily extracted from section 2.2. For instance, for the repeated case at Ω₀ = π there are three linearly independent solutions with

$$\begin{equation} u_0=f^{(1)}({\bf X})\sin\pi \xi_1 +f^{(2)}({\bf X})\sin\pi \xi_2 +f^{(3)}({\bf X})[\cos\pi \xi_1 -\cos\pi \xi_2]. \end{equation} \tag{ 40 }$$

Following through the methodology one arrives at 2Ω⁴₁ = (2π)²(κ²₁ + κ²₂) which gives the two linear asymptotics, and similarly at the other points. Table 1 summarizes the coefficients for the first two Dirac-like cones at point Γ. The asymptotics are shown in figure 7 and demonstrate that the linear behaviour is perfectly captured.

Table 1. The coefficients α_ii, necessary in equation (26) for the two Dirac-like cones at point Γ of the Brillouin zone for a doubly periodic array of constrained points.

α11	α22	Ω0
2π2	2π2	π
4π2	4π2	$\sqrt {2}\pi$

Table 2. The first six standing wave frequencies for in-phase waves at Γ, cf figure 4, together with associated values for T₁₁ and T₂₂. Symmetry between T₁₁ and T₂₂ is broken when the multiplicity of the eigenvalue is greater than two. Negative group velocity is demonstrated by the negative sign of both T_ij's.

T11	T22	Ω0
0.698 841 854 976 085	0.698 841 854 976 085	1.700 916 617 386 99
7.867 675 441 589 871	7.867 675 441 589 871	3.361 627 184 739 501
0.3230	−8.998	3.632 730 109 763 024
4.9681	14.2899	3.632 730 129 146 08
4.775 527 931 399 718	4.775 527 931 399 718	4.148 661 549 527 329
−4.279 843 054 160 115	−4.279 843 054 160 115	4.877 003 447 953 185

Table 3. The first five frequencies for in-phase standing waves and their respective T_ij coefficients for a doubly periodic array of constrained points (see figure 6). Equal T₁₁ and T₂₂ coefficients are for the single multiplicity eigenmodes. The second and fourth set of T_ij are for the quadratic modes of the Dirac-like cones.

T11	T22	Ω0
0.979 719 694 474 189	0.979 724 590 049 377	0.624 563 144 241 832
1	−22.34	π
12.439 613 406 598 808	12.439 765 169 617 527	3.388 815 608 872 719
−16.77	18.77	$\sqrt {2}\pi$
18.856 092 144 598 160	18.855 405 067 087 230	4.670 779 386 112 907

4.4. Square inclusions

Our methodology is not limited to circular inclusions and is easily applied to any shape, we briefly consider square inclusions taking a square of sidelength $\sqrt {2}$ and see in figure 8 that the dispersion curves are not unlike those of the large cylinder in figure 4. An isolated lowest mode separating a zero-frequency stop band from a wide stop band, yet another stop band arises near Ω = 6. The asymptotics from HFH again reproduce the behaviour near the edges of the Brillouin zone and can be used to construct effective HFH equations. The HFH results in figure 4 are almost completely indistinguishable throughout the lowest two branches and are highly accurate near the edges of the Brillouin zone for the other branches. Interestingly, the orientation of the square matters strongly and rotating it progressively flattens the lowest (lower frequency) dispersion curves until ultimately, for a rotation of π/4, they become completely flat leading to resonant states. The rotation gradually isolates each piece of material from its neighbours leaving resonant blocks that vibrate with the eigenfrequency of Dirichlet squares of sidelength $\sqrt 2$ $\Omega ^2=(m\pi /\sqrt 2)^2+(n\pi / \sqrt 2)^2$ with m,n non-zero integers.

Zoom In Zoom Out Reset image size

4.4.1. Standing wave modes

It is noticeable in figures 8 and 9(a) that completely flat modes exist which represent standing waves for a whole range of Bloch wavenumbers. The eigenfunctions, for the non-tilted case, for the first four standing wave frequencies, are shown in figure 10; it is clear that the field is concentrated along parallel and opposite sides and this suggests that a simple equivalent waveguide model can be constructed. Taking a rectangle with a Neumann condition on one side and Dirichlet conditions on the other three sides where the Neumann condition on one side is necessary to cover the appropriate symmetry. The length and width of the rectangle depend on the tilt of the square hole. For the non-tilted case the length is taken to be the length of a cell and the width as the half width of the cell minus the inner square. The resulting frequency estimate Ω_E reads

$$\begin{equation} \Omega_E=\pi\left(\left(\frac{2n-1}{l-w}\right)^2+\left(\frac{m}{l}\right)^2\right)^{\frac{1}{2}} \quad{\rm for}\quad n,m \in \mathbb N^*, \end{equation} \tag{ 41 }$$

where l and w/2 are, respectively, the effective length and width of the waveguide. For the non-tilted square hole l₀ = 2 and $w_0=2-\sqrt 2$ where for a tilt of π/8 they reduce to the values of l_π/8 = 1.6051 and w_π/8 = 0.4335. Equation (41) then leads to the 'x' crosses placed in figures 8 and 9(a) which are good predictions of the standing wave frequencies, and the approximate eigensolutions are shown in figure 10. The precision of the estimates degenerates as the frequency increases as the field near the ends of the imaginary rectangle gradually leaks into the undisturbed ligaments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now connect the HFH theory with the exciting phenomena that are topical in photonics. Once again, the matrix properties of the photonic crystals are a = ρ = 1.

5.1. Cloaking effects near Dirac-like cones

As noted in recent papers (Chan et al 2012) a curious phenomenon occurs in photonics near Dirac-like points which are the triple crossings shown in, for instance, figure 7. Let us consider the lowest frequency triple point, near wavenumber M, shown in figure 6; for the full finite element numerical simulations we actually consider very small cylinders of radius 0.01 in a square array. The triple crossings occur because the dispersion relation of a perfect medium (no cylinders) consists of light lines folded in space (one folded light line is shown in figure 6). The folded light lines happen to satisfy the boundary condition for the infinitesimal point holes, and are perturbed slightly for finite cylinders (as in figure 7); the middle curve passing between the Dirac-like cone is created by the inclusion and so these effects coexist. Importantly, this means that there is near perfect transmission through an array of cylinders close to these triple crossing frequencies as shown in figures 11(a) and (b) as the incoming plane wave, at Ω = 2.08, does not see the inclusions. Inserting a large Dirichlet hole within the array leads to virtually no reflection, the plane wave is hardly reflected, and its transmission is nearly perfect (figure 11(c)) except for a slight shadow zone in clear contrast to the uncloaked Dirichlet hole (figure 11(d)). If we remove the array of small holes, the defect strongly scatters the wave.

Zoom In Zoom Out Reset image size

Replacing the array of holes with HFH is shown in figure 12 wherein the quasi-periodic medium is replaced by a homogeneous medium with the effective properties given by HFH. The excitation is at a frequency Ω = 2.08, close to the standing wave frequency Ω₀ = 2.22. Inserting the values for l = 1, Ω₀, Ω and β = 2/π² into equation (32) we obtain the effective continuum equation

$$\begin{equation} {\skew5\hat{f}}_{0,x_ix_i}({\bf x})+0.0764{\skew5\hat{f}}_{\!0}({\bf x})=0, \end{equation} \tag{ 42 }$$

where an effective refractive index can be extracted as explained in section 3 to yield η_eff = 0.0177. Simulations show qualitatively the same effects, one has transmission through the slab with plane waves incoming.

Zoom In Zoom Out Reset image size

5.2. Lensing and wave guiding effects

Returning to the introduction we show in figures 2(a) and (c) a PC composed of an array of 196 holes with radius 0.4 and a centre-to-centre spacing of 2; the corresponding dispersion curves are shown in figure 4.

We now interpret these strongly anisotropic beams and how the HFH represents this. In figures 2(a) and (b) we choose frequency Ω = 1.98, which from figure 4 is near the standing wave frequency Ω₀ = 1.966 at point X(π/2,0). The reciprocal space is only shown for a portion of the irreducible Brillouin zone and one can use reflections of the triangle chosen to fill the entire square in the Brillouin zone; due to these internal symmetries the effect is not only due to the first mode of figure 4 at point X(π/2,0) of the reciprocal space but also to the first mode, at point G(0,π.2) (not shown). The effect is reproduced in figure 2(b) by combining two effective medium equations, one for point X and one for point G, each one responsible for a single diagonal, and the HFH equations are

$$\begin{equation} -1.4778{\skew5\hat{f}}_{0,x_1x_1}({\bf x})+0.8837{\skew5\hat{f}}_{0,x_2x_2}({\bf x})+(\Omega^2-\Omega_0^2){\skew5\hat{f}}_{\!0}({\bf x})=0, \end{equation} \tag{ 43 }$$

$$\begin{equation} 0.8837{\skew5\hat{f}}_{0,x_1x_1}({\bf x})-1.4778{\skew5\hat{f}}_{0,x_2x_2}({\bf x})+(\Omega^2-\Omega_0^2){\skew5\hat{f}}_{\!0}({\bf x})=0 \end{equation} \tag{ 44 }$$

with (Ω² − Ω²₀) = 0.0552. Note that these effective equations have opposite signs in the coefficient, and so the governing equations are hyperbolic and not elliptic; the result is that these direct energy along rays or characteristics and the angle is given by the relative magnitude of T₁₁ and T₂₂; in this case roughly equal with the energy directly along the diagonals.

To illustrate this further in panel (c) of figure 2 we show a Saint George's cross effect obtained near the standing wave frequency Ω₀ = 2.744. By homogenizing the PC and exciting it at a frequency of Ω = 2.75 we obtain two effective medium equations, one for X and one for G which are

$$\begin{equation} 3.1094{\skew5\hat{f}}_{0,x_1x_1}({\bf x})+0.085{\skew5\hat{f}}_{\!0,x_2x_2}({\bf x})+(\Omega^2-\Omega_0^2){\skew5\hat{f}}_{\!0}({\bf x})=0, \end{equation} \tag{ 45 }$$

$$\begin{equation} 0.085{\skew5\hat{f}}_{0,x_1x_1}({\bf x})+3.1094{\skew5\hat{f}}_{\!0,x_2x_2}({\bf x})+(\Omega^2-\Omega_0^2){\skew5\hat{f}}_{\!0}({\bf x})=0. \end{equation} \tag{ 46 }$$

Each equation is responsible for a single line of the cross. The orientation of this cross is explained by the almost zero coefficient, in each effective equation, that does not permit propagation in the vertical and horizontal directions for the respective equations (45) and (46). The coefficient in front of ${\skew5\hat{f}}_{\!0}({\bf{ x}})$ is equal to 0.033. It is clear therefore that the strong anisotropy that is witnessed in full FE numerical simulations has its origins in the size and sign of the T_ij coefficients from HFH and the effective equations can be used to gain quantitative understanding and predictions.

Another striking application involving wave propagation and an array of small holes of radius 0.01 is shown in figures 13–15 wherein we have tilted the array through an angle π/4. The dispersion curves in figure 6 are for holes of radius exactly zero, but provide very good comparison for this small holes case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We obtain a highly directed emission at frequency Ω = 1.6 near the standing wave frequency Ω₀ = π/2, shown in figure 14(a) with its effective medium counterpart in figure 14(b); in the dispersion curves this is the lowest standing wave frequency at X. There are two effective medium equations at this frequency which yield

$$\begin{equation} -11.42{\skew5\hat{f}}_{0,x_1x_1}({\bf x})+{\skew5\hat{f}}_{\!0,x_2x_2}({\bf x})+(\Omega^2-\Omega_0^2){\skew5\hat{f}}_{\!0}({\bf x})=0, \end{equation} \tag{ 47 }$$

$$\begin{equation} {\skew5\hat{f}}_{0,x_1x_1}({\bf x})-11.42{\skew5\hat{f}}_{\! 0,x_2x_2}({\bf x})+(\Omega^2-\Omega_0^2){\skew5\hat{f}}_{\!0}({\bf x})=0, \end{equation} \tag{ 48 }$$

where the difference of the squared frequencies is equal to 0.0926. For forcing within the slab one again has a strongly anisotropic response, the waves within the medium form a cross-shape which, when it strikes the edge of slab, all gets radiated out into the surrounding medium. Placing a source outside the slab, again at a frequency near Ω₀ = π/2, figure 15(a) shows a PC that behaves like a flat lens. In panel (b) of that figure the effective medium yields the same effect and is also governed by two equations of the form of (30) with T_ij coefficients equal to those of equation (48). The frequency at which the lensing effect is obtained is Ω = 1.56 and so the coefficients in front of ${\skew5\hat{f}}_{\!0}({\bf{ x}})$ are equal to −0.0338 in this case.

Figures 13(a) and 14(c) show lensing effects where the PC is excited at Ω = 2.7, this frequency is chosen as follows. When the array is rotated the relevant light line in the surrounding material, and its folding, emerge from the point M and are shown in figure 6; we see that the energy from the source couples (as the light line crosses the full dispersion curves at Ω = 2.7) into one of the linear Dirac-like cone modes that passes through Ω₀ = π; this is the standing wave frequency of the Dirac-like cone at Γ. We therefore model the effective medium with HFH using the analysis near Ω = π, and it yields very similar results in figures 13(b) and 14(d), by homogenizing the PC with an effective material governed by the equation,

$$\begin{equation} {\skew5\hat{f}}_{0,x_ix_i}({\bf x})+0.3371{\skew5\hat{f}}_{\!0}({\bf x})=0 \end{equation} \tag{ 49 }$$

obtained by simply replacing the source and standing wave frequencies together with β = 1/(2π²) in equation (32).

Finally, figure 16(a) shows a unidirective PC for the first two flat modes of the band diagram in figure 9(a). The standing wave frequencies of the two flat modes in question are Ω₀ = 7.6 and Ω₀ = 8.28. The sources are excited near these frequencies at Ω = 7.7 and Ω = 8.25. HFH yields effective medium equations for each of the modes that are, respectively,

$$\begin{equation} \quad9.3649{\skew5\hat{f}}_{0,x_1x_1}({\bf x})+(\Omega^2-\Omega_0^2){\skew5\hat{f}}_{\!0}({\bf x})=0 \end{equation} \tag{ 50 }$$

and

$$\begin{equation} \quad-45.8434{\skew5\hat{f}}_{0,x_1x_1}({\bf x})+(\Omega^2-\Omega_0^2){\skew5\hat{f}}_{\!0}({\bf x})=0. \end{equation} \tag{ 51 }$$

Equations (50) and (51) clearly show the unidirectivity of the two effective media. Similar equations with interchanged coefficients stand for the symmetric location of the Brillouin zone G(0,π/2). Panels (b)–(d) contain PCs with rotated square holes by an angle of π/4 with an unusual band diagram seen in figure 9(b). Panel (b) shows the guiding of a wave by means of perfect reflection for a line source excited just outside the standing wave frequency of Ω = 4.95 while in panel (c) the same PC but this time rotated by an angle of π/2 acts as a perfect reflector. Panel (d) finally excites the medium right on the standing wave frequency of Ω = 5 and standing waves appear between the holes. Small gaps between the hole's edges make it possible to channel energy from one cell to the other.

Zoom In Zoom Out Reset image size

In this paper, we have presented a range of applications of the HFH theory in the context of zero-frequency stop band structures, which were previously thought to be non-homogenizable: HFH has succeeded in capturing the finer details of both dilute and densely packed photonic and phononic crystals. Striking physical effects such as cloaking via Dirac-like cones and directive antennae and endoscope effects have been provided, and fully explained, via HFH; importantly this is all shown to be related to the generally anisotropic T_ij coefficients that encode the local behaviour of the effective medium or their equivalent terms for multiple crossings. The range of unsolved homogenization problems in the wave community is vast, and this new method now opens the door to efficient simulation of multiscale structures.

The current theory is limited to periodic, or nearly periodic media; however, one can envisage extensions in a variety of different directions. In this paper we deal with photonic and phononic periodic structures, but it should be possible to adapt our results to quasi-crystals using a generalized form of the Floquet–Bloch theorem in upper-dimensional spaces. Stochastic cases where the hole positions involve some random perturbation, say, might prove more arduous from a mathematical standpoint although a physicist might argue the period would be replaced by the mean distance between the inclusions. The key point in order to relax the constraint of classical homogenization is the assumption that the Floquet–Bloch theorem is applicable at least to leading order. Whenever this can be done, or an extension of the Floquet–Bloch theorem can be used, HFH is applicable. Also of interest are diffusion problems in composites, classical homogenization is here again constrained by low frequencies (e.g. of a periodic heat source), whereas much of the exciting physics lies beyond the quasi-static regime again we hope that the insight provided by HFH will lead to progress in these directions.

The authors thank Julius Kaplunov, Kirill Cherednichenko, Aleksey Pichugin and Evgeniya Nolde for many enjoyable and useful conversations. RVC thanks the EPSRC (UK) for support through research grant number EP/J009636/1 and the University of Aix-Marseille for support via a Visiting Professorship. SG is thankful for an ERC starting grant (ANAMORPHISM) that facilitates collaboration with Imperial College London.