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

We present an accurate methodology for representing the physics of waves, for periodic structures, through effective properties for a replacement bulk medium: This is valid even for media with zero frequency stop-bands and where high frequency phenomena dominate. Since the work of Lord Rayleigh in 1892, low frequency (or quasi-static) behaviour has been neatly encapsulated in effective anisotropic media. However such classical homogenization theories break down in the high-frequency or stop band regime. Higher frequency phenomena are of significant importance in photonics (transverse magnetic waves propagating in infinite conducting parallel fibers), phononics (anti-plane shear waves propagating in isotropic elastic materials with inclusions), and platonics (flexural waves propagating in thin-elastic plates with holes). Fortunately, the recently proposed high-frequency homogenization (HFH) theory is only constrained by the knowledge of standing waves in order to asymptotically reconstruct dispersion curves and associated Floquet-Bloch eigenfields: It is capable of accurately representing zero-frequency stop band structures. The homogenized equations are partial differential equations with a dispersive anisotropic homogenized tensor that characterizes the effective medium. We apply HFH to metamaterials, exploiting the subtle features of Bloch dispersion curves such as Dirac-like cones, as well as zero and negative group velocity near stop bands in order to achieve exciting physical phenomena such as cloaking, lensing and endoscope effects. These are simulated numerically using finite elements and compared to predictions from HFH. An extension of HFH to periodic supercells enabling complete reconstruction of dispersion curves through an unfolding technique is also introduced.


Introduction
Over the past 25 years, many significant advances have created a deep understanding of the optical properties of photonic crystals (PCs) (Yablonovitch 1987, John 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 (Pendry et al. 1996, Nicorovici et al. 1995. 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 directions. This effect is well known (Joannopoulos et al. 2008, Zolla et al. 2005 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 in regard with its incidence (Markos & 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. It is therefore possible to create seismic shields for low-frequency waves (Brule et al. 2013, Antonakakis et al. 2013b. Extending these ideas to higher dimensions allows for a much greater range of optical effects: all-angle negative refraction (Luo et al. 2002, Zengerle 1987, ultrarefraction (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 medium can be created that reproduces these effects.
The 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) for cylinders, to name but a few, that allow one to visualize the various effects and design PCs. The multipole method takes its root in a seminal paper of Lord Rayleigh (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, Figure 1. High Frequency Homogenization (HFH) principle: An elementary cell (a) of sidelength 2l, modelled by a fast oscillating variable ξ, is repeated periodically within a supercell (b) of sidelength L, modelled by a slow variable X, which is itself repeated periodically in space (c). One then assumes that the parameter = l/L is small, and its vanishing limit thereafter studied. The leading order, homogenized, term of Floquet-Bloch eigenfields within the crystal are then sought as u 0 (X, ξ) = f 0 (X)U 0 (ξ; Ω 0 ), wherein f 0 accounts for variations of the fields on the order of the supercells, and U 0 captures their fast oscillations in the much smaller cells, when either periodic or antiperiodic conditions are enforced on the cells: Perturbing away from these standing waves of frequency Ω 0 allows for a complete reconstruction of the Bloch spectrum and associated Floquet-Bloch eigenfields.
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 a linear behaviour in the neighbourhood of zero frequency, except of course in the singular case whereby Dirichlet data is 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 do 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 past 35 years by applied mathematicians primarily for solving partial differential equations, with rapidly oscillating periodic coefficients, in the context of thermostatics, continuum mechanics or electrostatics (Bensoussan et al. 1978, 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 & Zolla 2000, Wellander & Kristensson 2003 and high-contrast two-dimensional (Zhikov 2000, Bouchitté & Felbacq 2004, Cherednichenko et al. 2006) photonic crystals, that have important potential applications in photonics. 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 & 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 on the order to −2 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 (O'Brien & 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 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 homogenisation (Conca et al. 1995, Allaire & Piatnitski 2005, Birman & Suslina 2006, Hoefer & 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 homogenisation limitations.
In this paper, we apply the HFH theory proposed by one of us three 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 two years ago . 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 & Nicorovici 2006, Guenneau et al. 2007) to anti-earthquake seismic shields (Antonakakis et al. 2013a, Brule et al. 2013. We show in Fig. 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's 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 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 Fig. 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   , this strong frequency dependent anisotropy is a recurrent feature of PCs. One sees that the HFH neatly reproduces in Fig. 2(b,d) the fine features of the full finite element simulations. Interpretation is provided by using the Brillouin zone of Fig. 3 and dispersion curves shown in Fig. 4. The square array of holes behaves as two dramatically contrasting effective media at frequency Ω = 1.966, see Fig. 2(a,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, wherein the dispersion curve is nearly flat is the M X direction; both these frequencies are well beyond the range of applicability of classical homogenization.
Our aim is to generate HFH for arrays of holes with Dirichlet conditions, so this is for TM waves in electromagnetism, and then compare through full numerical simulations of HFH and finite elements to show how the behaviour of the numerous interesting optical effects emerge naturally through the coefficients of our effective medium. This is intertwined with a 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 is key. In section 2 we set up the mathematical framework leading to the effective medium equation (15), this is in a non-dimensional 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 these 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, and we use both circular and square holes for illustration, as shown in section 4. Local to 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 apriori, 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.

General theory
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 Fig. 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 article our emphasis is rather different, and is upon TM polarization for which the Dirichlet conditions induce different effects, such as zerofrequency stop bands, and require a modified theory.
Assuming an exp(−iωt) time dependence that is considered understood, and henceforth suppressed, the governing equation is, 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 1 /l, x 2 /l), namelyâ(ξ) andρ(ξ). In the context of photonics, these are the inverse of permeability (µ −1 ), and permittivity (ε), respectively, which are related to the wavespeed of light in vacuum c via the refractive index n as follows: εµ = n 2 c −2 . 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ω −1 µ −1 ∇×E. In the context of phononics,â(ξ) andρ(ξ) 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 & Guenneau 2013) metamaterials. The length of the direct lattice base vectors are taken equal and to be 2l, and define a short lengthscale. 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,â ≡â 0 a(ξ) and ρ ≡ρ 0 ρ(ξ) is convenient. Equation (1) then becomes as the non-dimensional frequency andĉ 0 = â 0 /ρ 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, For wave propagation through a perfect lattice there exist standing waves at specific eigenfrequencies, Ω 0 , and their eigenmodes satisfy periodic, or anti-periodic (out-ofphase) 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 where u ,ξ i denotes partial differentiation with respect to ξ 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 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 We start by expressing a solution for the leading order equation. For a specific eigenvalue Ω 0 there is a corresponding eigenmode U 0 (ξ; Ω 0 ), periodic in ξ, and the leading order solution is The function f 0 (X) is determined later on, indeed the whole point is to deduce a longscale continuum partial differential equation for f 0 entirely upon the long-scale. We begin with a treatment of isolated eigenvalues, however repeated eigenvalues also arise, and require 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 2 , where S denotes the surface of the cell in the ξ coordinates, yield and are set in the short-scale ξ so, for i = 0, U 0 (ξ; Ω 0 )| ∂S 2 = 0.
To deduce the long-scale PDE for f 0 (X) we follow (Craster et al. 2010), making changes where necessary due to the change in boundary conditions. Equation (7) is multiplied by U 0 and integrated over the cell's surface to obtain The first integral of the right hand side is converted to a path integral along ∂S using a corollary of the divergence theorem. The periodic conditions on ∂S 1 and the homogeneous Dirichlet conditions of U 0 on ∂S 2 make it vanish. We continue by subtracting the integral over the cell of the product of equation (6) with u 1 /f 0 to obtain Using Green's theorem equation (12) becomes where ∂S = ∂S 1 + ∂S 2 . Due to periodic (or anti-periodic) boundary conditions of u and its first derivatives with respect to ξ i on ∂S 1 , and due to the homogeneous Dirichlet boundary conditions of the u i 's on ∂S 2 , the left hand side of equation (13) vanishes. It follows that Ω 1 = 0, which is an important deduction as it implies that the asymptotic behaviour of dispersion curves near isolated eigenfrequencies is at least quadratic. We then solve for u 1 (X, ξ) and obtain The first term is simply a homogeneous solution that is irrelevant, and the vector function U 1 is a particular solution. Quite often exact solutions for u 1 or U 1 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 0 (X) and Ω 2 2 . Similarly to the first order equation, we multiply equation (8) by U 0 then subtract the product of equation (7) with u 2 /f 0 , finally we integrate over the cell's surface to obtain the following effective equation for f 0 (15) which is the main result of this article. In particular, the coefficients T ij encode the anisotropy at a specific frequency and the t ij 's are integrals over the small-scale cell where U 1 i is the ith component of vector function U 1 . Note that there is no summation over repeated indexes for t ii . The PDE for f 0 , equation (15), is crucial as the local microstructure is completely encapsulated within the tensor T ij ; these are, for a specific structure at an Ω 0 , 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 or when specific effects occur. The structure of the tensor depends upon the boundary conditions of the holes, the results here, for instance, are 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 (15) is simple and quick to solve numerically, or even by hand.
One primary aim of the present article 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 (X) = exp(iκ j X j / ). In this notation κ j = K j − d j and d j = 0, π/2, −π/2 depending on the location in the Brillouin zone. Equation (15) gives and these locally quadratic dispersion curves are completely described by Ω 0 and the tensor T ij .

Classical singularly perturbed zero-frequency limit
It is natural, at this point, to contrast with usual homogenization theories. As is wellknown (Nicorovici et al. 1995, McPhedran & Nicorovici 1997 one cannot homogenize the TM polarized case because the conventional approach only works at low frequencies in a quasi-static limit. Some authors attribute this to a singular perturbation, or noncommuting limit, problem whereby letting first the Bloch wavenumber and then the frequency tend to zero, or doing it the way around, leads to a different result (Poulton et al. 2001). The approach used in this article overcomes this by releasing the theory from the low-frequency constraint. For TE polarization (Antonakakis et al. 2013a) one can explicitly connect the theories and show that the low frequency theory is a subset of the high frequency approach. The TM case differs fundamentally as there is a zero-frequency stop-band. We now prove that the usual homogenization cannot give the asymptotic behaviour even at low-frequency. Setting Ω 2 = 2 Ω 2 2 + . . ., that is, going to low frequency, and assuming that u 0 (X, ξ) = f (X)U 0 (ξ) is non-zero one arrives at a contradiction: At leading order U 0 (ξ) satisfies ∇ 2 U 0 = 0 where for brevity we have set a = 1, complemented by U 0 = 0 on the hole and periodic (anti-periodic) boundary condition on the edge of the cell. From the uniqueness properties of the Laplacian U 0 ≡ 0 is the unique solution and hence usual homogenization cannot find the asymptotics at low frequency in this Dirichlet case.

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, where we sum over the repeated superscripts (l). Proceeding as before, we multiply equation (7) by U There is now a system of coupled partial differential equations for the f (l) 0 and, provided Ω 1 = 0, the leading order behaviour of the dispersion curves near the Ω 0 is now linear (these then form Dirac-like cones). These coupled partial differential equations on the long-scale now replace (15) near these frequencies.
For the perfect lattice and Bloch wave problem, we set f 0 exp(iκ j X j / ) and obtain the following equations, where, and The system of equations (21) is written simply as, with C ll = Ω 2 1 B ll and C ml = iκ j A jml / for l = m. One must then solve for Ω 2 1 = ± √ α ij κ i κ j / when the determinant of C vanishes and the asymptotic relation is, If Ω 1 is zero, one must go to the next order and a slightly different analysis ensues.

Repeated eigenvalues: Quadratic asymptotics
Assuming that Ω 1 is zero, 1 k (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) 0 and u 2 ((aU 0,ξ i ) ,ξ i + Ω 2 0 ρU 0 ) and then integrating over the elementary cell gives as a system of coupled PDEs. The above equation is presented more neatly as For the Bloch wave setting, using f (l) 0 (X) =f (l) 0 exp(iκ j X j / ) we obtain the following system, .., p (28) and this determines the asymptotic dispersion curves.

Effective dispersive media
The major goal of HFH is to represent the periodic medium of finite extent using an effective homogeneous medium and the main result of this article is in achieving this, and for clarity this is summarized in this section back in the dimensional setting.
Transforming equation (15) back to the original x i = X i L coordinates and using the solution of Ω 2 we obtain an effective medium equation forf 0 that is, The nature of equation (29) is not necessarily elliptic since T 11 and T 22 can take different values and/or different signs. In the illustrations herein the T ij = 0 for i = j. The hyperbolic behaviour of equation (29) 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 the 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.
At Dirac-like cones the linear behaviour of the effective medium yields an equation slightly different from (29). Regarding the zero radius holes, presented in the following section 4.3, it consists of a system of three coupled PDE's that uncouple to yield one identical PDE for all f (i) 0 s that is of the form, where β is a coefficient equivalent of the T ij but this time is the same for all combinations  Table 2. The first six standing wave frequencies for in-phase waves at Γ, cf. Fig. 4, together with associated values for T 11 and T 22 . Symmetry between T 11 and T 22 is breaking when the multiplicity of the eigenvalue is greater than two. Negative group velocity is demonstrated by the negative sign of both T ij 's.
of i and j. The quadratic mode that emerges in the middle of the linear ones follows equation (29).

Dispersion curves
We now generate dispersion curves for circular and square holes 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.

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 (Zolla et al. 2005, Movchan et al. 2002.

Large circular holes
The full dispersion curves are computed numerically using COMSOL (finite element code (COMSOL 2012)) and then asymptotically with the HFH using the equations developed in section 2. Fig. 4 illustrates the typical dispersion curves versus the asymptotics for wavenumbers, on the specified path of Fig. 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 11 , T 22 values for point Γ of the Brillouin zone. Note how T 11 = T 22 for all single eigenfrequencies, but that symmetry breaks for the double roots. Moreover the signs of the T 11 and T 22 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.

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 asymptotically coverage, see Fig. 5. As an example, we take an elementary square cell of length 2 with circular holes of radius 0.4, section 4.1 and  Table 3. The first five frequencies for in-phase standing waves and their respective T ij coefficients for a doubly periodic array of square cells with zero radius holes (see Fig. 6). Equal T 11 and T 22 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. 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 ξ 1 direction with anti-periodic boundary conditions on ξ 1 but periodic on ξ 2 and finally a rectangular macro cell composed of two elementary cells in the ξ 2 direction with anti-periodic conditions in ξ 1 and ξ 2 . The asymptotics at points M , I, J as well as the usual points Γ, M and X are illustrated in Fig. 5.

Circular inclusions of small radius
As the hole radius decreases, the dispersion curves gradually transition toward those of the infinitesimal holes shown in Fig. 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.
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 toward those points to replace the quadratic behaviour of all the multiple modes, except one, as seen in Fig. 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 hep-tuple crossing. Reassuringly, the asymptotic theory captures these multiple modes whether the crossings are double, triple, or more as is illustrated in Figs. 4, 6 where the circular holes are respectively of radius 0 and 0.4. For the multiple crossing points with multiplicity higher than two, both linear and quadratic dispersion curves arise. One can proceed, as in section 2.2, to obtain the different Ω 1 coefficients and one of them is zero; proceeding to higher order and apply equations (26) to (28) gives the quadratic behaviour of the middle mode emerging from the triple crossings. Again the T ii contain the 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 zero radius and for 0.1 radius (not shown) geometries are almost indistinguishable by 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 almost touching standing wave frequencies depends on the size of the inclusion and tends to zero as the inclusion shrinks down to a point: Fig. 7 illustrates a close-up of a triple point comparing inclusions of radius 0 and 0.1. This is actually important for the Dirac-like cones and associated effects.

Zero radius, Fourier series and Dirac-like cones
A doubly periodic array of constrained points (circles of zero radius) is an attractive limit and has an immediate solution in terms of Fourier series as with κ = (κ 1 , κ 2 ) and N = (n 1 , n 2 ) and Ω and κ are related via the dispersion relation This is valid provided the double sum does not have a singularity, which would occur whenever 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 Fig. 6. It is also clear from the dispersion relation that Ω = 0 is not a solution at κ = 0 and hence that there is a zero-frequency stopband and conventional long-wave homogenization is doomed to failure cf. section 2.1. A set of dispersion curves are shown in Fig. 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 are called generalized Dirac-like points: Generalized Dirac-like points 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 Ω 2 = π 2 m such that m = n 2 1 + n 2 2 with m integer; the multiplicity depends on the ways in which n 1 and n 2 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 2 + 5 2 = 7 2 + 1 2 or the sum of two squares can also be a perfect square, i.e., m = 25 = 5 2 = 4 2 + 3 2 . Equation u ,x i x i + π 2 (n 2 1 + n 2 2 )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 i x + n j y)), sin(π(n i x − n j y)), cos(π(n i x + n j y))-cos(π(n p x − n k y)) 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 Ω 2 = π 2 with n i = 1, n j = 0 has multiplicity 3 where Ω 2 = 5π 2 has multiplicity 7 and Ω 2 = 50π 2 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 1 ) 2 + n 2 2 , it is possible to obtain solutions of multiplicity one as in the case of Ω 0 = π/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 0 from (31) and deduce U 1 = (U 11 , U 12 ) as and thus one can extract the T ij in (15) by doing the integrals by hand, the off-diagonal terms are zero, and T 11 and T 22 are and an identical equation for T 22 but with 1 and 2 interchanged, where (κ 1 , κ 2 ) 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 Fig. 6. The linear cases at the triple crossings are easily extracted from section 2.2. For instance, for the repeated case at Ω 0 = π there are three linearly independent solutions with Following through the methodology one arrives at 2Ω 4 1 = (2π) 2 (κ 2 1 + κ 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 Fig. 7 showing that the linear behaviour is perfectly captured.

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 side √ 2 and see in Fig. 8  Figure 8. The dispersion curves for square inclusions, of side √ 2, from numerical simulation (solid) and from HFH theory (dashed). The HFH and numerics for the lowest two curves are virtually indistinguishable. The crosses on the frequency axis are frequencies predicted by solving Helmholtz's equation in a waveguide consisting of a rectangle with Neumann data on one side and Dirichlet data on the other sides, see also Figure. 10. Their values in increasing order are 5.59, 6.22, 7.14 and 8.26. that the dispersion curves are not unlike those of the large cylinder in Fig. 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 Fig. 4 are almost completely indistinguishable all along 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 side √ 2 Ω 2 = (mπ/ √ 2) 2 + (nπ/ √ 2) 2 with m, n non-zero integers.
4.4.1. Standing wave modes It is noticable in Figs. 8 and 9a 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 Fig. 10; it is clear that the field is concentrated along parallel and opposite sides and this suggests 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 depends 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, where l and w/2 are respectively the effective length and width of the waveguide. For the non-tilted square hole l 0 = 2 and w 0 = 2 − √ 2 where for a tilt of π/8 they reduce to the values of l π/8 = 1.6051 and w π/8 = 0.4335. Equation (37)   (c) Same as in panel (b) when a clamped rectangular obstacle is placed inside the array; (d) Plane wave incident from above on the clamped obstacle at frequency Ω 0 for comparison; The slightly reduced amplitude in forward scattering in (c), but the lack of backward scattering, is a hallmark of cloaking.
frequencies, and the approximate eigensolutions are shown in Fig. 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.

Lensing, guiding and cloaking effects
We now connect the HFH theory with the exciting phenomena that are topical in photonics. Figure 12. Panel (a) shows nearly perfect transmission for a plane wave at frequency Ω 0 incident from above on an effective medium with properties computed by HFH for a frequency near Ω 0 = 2.215 with PML on either sides of the computational domain; A Plane wave incident from above on the clamped obstacle at frequency Ω 0 for comparison. Figure 13. Lensing effects in a PC (tilted array) of small Dirichlet holes: (a) Lensing for a line source inside the PC at frequency Ω = 2.7; (b) Using equation (29) the PC is replaced by an effective medium in order to yield the same effect as in (a).

Cloaking effects near Dirac-like cones
As noted in recent articles (Chan et al. 2012) a curious phenomenon occurs in photonics near Dirac-like points which are the triple crossings shown in, for instance, Fig. 7. Let us consider the lowest frequency triple point, near wavenumber M , shown in Fig. 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 Fig. 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 Fig.  7); the middle curve passing between the Dirac-like cone is created by the inclusion and so these effects co-exist. Importantly, this means that there is near perfect transmission The PC of (a) is replaced by an effective medium obtained by HFH with equation (29) in order to obtain the same physical effect; (c) Endoscope for a line source at Ω = 2.7; (d) The PC of (c) is replaced by an effective medium using equation (30) to obtain the same physical effect.
through an array of cylinders close to these triple crossing frequencies as shown in Fig.  11(a,b) as the incoming plane wave, at Ω = 2.215, 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 Fig. 11(c) except for a slight shadow zone in clear contrast to the uncloaked Dirichlet hole Fig. 11(d); If we remove the array of small holes, the defect is strongly scattered by the wave.
Replacing the array of holes with HFH is shown in Fig. 12 the quasiperiodic medium is replaced by a homogeneous medium with the effective properties given by HFH. The excitation is at a frequency Ω = 2.12, close to the standing wave frequency Ω 0 . Inserting the values for l = 1, Ω 0 , Ω and β = 2/π 2 into equation (30) we obtain the effective continuum equation, Simulations show qualitatively the same effects, one has transmission through the slab with plane waves incoming.

Lensing and wave guiding effects
Returning to the introduction we show in Fig. 2(a,c) a PC composed of an array of 196 holes with radius 0.4; the corresponding dispersion curves are shown in Fig. 4. We now interpret these strongly anisotropic beams and how the HFH represents this. In Fig 2 (a,b) we choose frequency Ω = 1.98, which from Fig. 4 is near the standing wave frequency Ω 0 = 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 Fig. 4 at point X(0, π/2) of the reciprocal space, but also to the first mode, at point G(π/2, 0) (not shown). The effect is reproduced in Fig. 2 (b) by combining two effective medium equations, one for point X and one for point G, each one is responsible for a single diagonal, and the HFH equations are with (Ω 2 − Ω 2 0 ) = 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 the T 11 and T 22 ; in this case roughly equal with the energy directly along the diagonals.
To illustrate this further in panel (c) of Fig. 2 we show a St George's cross effect obtained near the standing wave frequency Ω 0 = 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 that are, 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 direction for the respective equations (41) and (42). The coefficient in front off 0 (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 Figs. 13, 14 , 15 wherein we have tilted the array through an angle π/4. The dispersion curves in Fig. 6 are for holes of radius exactly zero, but provide very good comparison for this small holes case.
We obtain a highly-directed emission at frequency Ω = 1.6 near the standing wave frequency Ω 0 = π/2, shown in Fig. 14(a) with its effective medium counterpart in 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, where the difference of the frequency squares is equal to 0.0926. For forcing within the slab one again has a strongly anisotropic response, the waves within the medium forming 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 Ω 0 = π/2, Fig. 15(a) shows a PC that behaves like a flat lens. In panel (b) of that Fig. the effective medium yields the same effect and is also governed by two equations of the form of (29) with T ij coefficients equal to the ones of equation (44). The frequency at which the lensing effect is obtained is Ω = 1.56 and so the coefficients in front off 0 (x) are equal to −0.0338 in this case. Figs. 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 Fig. 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 Ω 0 = π; 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 Figs. 13(b) and 14(d), by homogenizing the PC with an effective material governed by equation, Figure 16. Lensing and guiding effect through a PC composed of square Dirichlet holes which we vary their tilt angle: (a) Lensing effects in a PC composed of Dirichlet square holes tilted by π/8 near the respective frequencies of 7.7 and 8.25 for the top and bottom figures (Fig 9a). (b) The dispersion curves in Fig. 9b show the band structure is composed of many band gaps separated by standing wave modes, so that the structure acts as a perfect reflector for any frequency in order to guide light. (c) shows a perfect reflector for a source exciting the medium at frequency 4.95, located in the second band gap of Fig. 8. (d) shows the same effect as in (c) for a frequency of 5.0 although in the present the filled parts of material are excited but no energy escapes the PC and acts as an energy trap.
obtained by simply replacing the source and standing wave frequencies together with β = 1/(2π 2 ) in equation (30). Finally, Fig. 16(a) shows a unidirective PC for the first two flat modes of the band diagram in Fig. 9(a). The standing wave frequencies of the two flat modes in question are Ω 0 = 7.6 and Ω 0 = 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, 9.3649f 0,x 1 x 1 (x) + (Ω 2 − Ω 2 0 )f 0 (x) = 0, and −45.8434f 0,x 1 x 1 (x) + (Ω 2 − Ω 2 0 )f 0 (x) = 0, Equations (46) and (47) clearly show the unidirectivity of the two effective media. Similar equations with interchanged coeffcients stand, for the symmetric location of the Brillouin zone G(0, π/2). Panels (b), (c) and (d) contain PCs with rotated square holes by an angle of π/4 with an unsual band diagram seen in Fig. 9(b). Panel (b) shows the guiding of a wave by means of perfect reflection for a line source excitated 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 for energy to pass from one cell to the other.

Concluding remarks
In this paper, we have presented a range of applications of the high frequency homogenization 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 antennas and endoscope effects have been provided, and fully explained, via HFH; importantly this is all shown to be related to the potentially anisotropic T ij coefficients that encode the local behaviour 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, classic 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.