## Abstract

Reflection and transmission of electromagnetic waves at the boundaries of periodic composites (electromagnetic/optical metamaterials) depends in general on both bulk and surface waves. We investigate the interplay of these two contributions using three-dimensional full-wave numerical simulations and a recently developed non-asymptotic homogenization theory.

© 2013 OSA

## 1. Introduction

The problem of homogenization of periodic electromagnetic composites has proved to be quite complicated and is currently the subject of lively research [1–7]. In many applications and simulations reported in the literature, the composite has the overall shape of a slab with elementary cells arranged on a cubic lattice with a period *h*. Let us assume that a slab of sufficiently large width *L* (that is, *L* = *Nh* where *N* ≫ 1) can be accurately described by some effective medium parameters. On the other hand, for a slab containing only one layer of elementary cells (*N* = 1), this description is expected to be substantially different or inapplicable. The composites used in practice will, most likely, have some intermediate value of *N*. This raises the questions of how large the number of elementary layers should be for the effective medium description to be applicable, what errors are incurred by using this description for a given finite *N*, and whether the effective medium parameters can depend on *N*. It should be noted that many simulations have revealed that the width effects are unimportant in slabs containing only a few layers of elementary cells [8–10]. These findings, however, are largely heuristic. It is desirable, therefore, to gain some qualitative understanding of the interplay between the sample width and the accuracy of the effective medium description.

The problem has been considered previously in a number of publications. In Ref. [11] it was shown that material parameters retrieved from the reflection and transmission coefficients of a composite slab (by means of the so-called S-matrix retrieval method [12, 13]) can, in fact, depend appreciably on the number of elementary cell layers used. In particular, parameters obtained for a single layer are inapplicable to a multilayer structure. Moreover, these parameters can depend on the angle of incidence of the illuminating plane wave and therefore cannot always be used in the traditional sense as characteristics of the material itself independently of the illumination conditions (see also relevant work in Ref. [8, 9, 14]).

In this paper, instead of using the transmission and reflection data for parameter retrieval, we investigate the problem by looking at the electromagnetic fields excited inside the composite. To this end, we use the homogenization methodology of Refs. [15–17]. In these references, the effective parameters are derived from the electromagnetic fields computed inside a given cell. In most practical applications, this computation is done numerically because analytical solutions are known only for a few rare special cases (e.g., one-dimensional periodically-layered medium). The advantage of this homogenization approach is that it allows one to study the effective parameters not only as functions of the number of layers, but also as functions of position relative to the slab surfaces. It is, in fact, the key objective of this paper to demonstrate that the effective parameters of surface layers can be defined independently from the bulk parameters and may differ from the latter. To keep numerical simulations and analysis relatively simple, we focus exclusively on the surface effects and do not consider a number of other relevant issues. In particular, even though the theory of Refs. [15, 17] applies to any illumination conditions, in the present paper we consider normal incidence only. Correspondingly, the effective parameters computed below represent the medium accurately for a normally-incident plane wave, but not necessarily for other types of illumination. The effective parameters we compute are guaranteed to apply to all types of illumination for sufficiently long waves, but this may cease to be true for shorter waves. We expect that the spatially-dispersive behavior (that is, inapplicability of a purely local effective medium theory) takes place in the spectral regions where the magnetic permeability is noticeably different from unity, e.g., for *λ*/*h* ≲ 5 [8, 18, 21] (*h* is the lattice period, see below). We note that the practical utility of nonlocal parameters, if these can be reasonably introduced at all, is debatable and, in any event, the onset of nonlocality tends to destroy all physical effects that enable the frequently discussed applications of metamaterials, such as super-resolution. We, therefore, do not investigate the short-wavelength spectral region or the effects of nonlocality in detail, but rather focus on the influence of the surface wave, which is present at both long and short waves.

We argue that the presence of a surface wave explains why the surface and bulk parameters are different and why the thickness of the sample affects the effective parameters obtained using the traditional transmission/reflection retrieval. The exponential decay of the surface wave away from the interface helps to explain why the dependence of parameters on the thickness of the slab is appreciable only for relatively thin samples, as observed in several previous studies [2, 4, 8, 9, 14]. We note that the surface wave exists due to the lack of discrete translational invariance in finite composites, and is conceptually similar to surface plasmon polaritons (SPPs), even though the typical physical realization of SPPs is somewhat different: the SPPs are conventionally excited on a metallic or polaritonic (e.g. SiC) surface separating two media with the opposite signs of the real part of *ε*.

A qualitative connection can be made between our study of surface parameters and transition layers whose importance has been emphasized by Simovski (e.g. [2]); he traces this subject back to the work of Drude. Our homogenization model allows one to rigorously define and quantify electromagnetic parameters of the transition layers.

The remainder of the paper is organized as follows. In Sec. 2 we briefly recount the theory of surface waves in finite, periodic, three-dimensional composites. In Sec. 3 we describe the homogenization approach used by us. In Sec. 4 we present the main results of this paper and, finally, Sec. 5 contains a brief discussion.

## 2. Surface waves – theory

In the case of smooth interfaces, SPPs can be excited only by an evanescent incident wave, i.e., from the near-field, or from a medium with a higher refractive index (in the Kretschmann geometry [19]). However, in periodically modulated media, SPPs can be exited by both propagating and evanescent incident waves and are, therefore, unavoidable. The excitation of SPPs in three-dimensional composites can be described in the framework of the generalized Ewald-Oseen theorem [20, 21]. The classical Ewald-Oseen theorem states that any wave incident on an interface with a homogeneous medium creates a wave of polarization; this, in turn, produces the extinction electromagnetic wave (which cancels exactly the incident wave) and a transmitted wave (which has exactly the same spatial dependence as the wave of polarization). The physical picture of refraction is thus made self-consistent from the microscopic point of view [21].

A similar situation exists in periodic composites but with an important variation. In this case, the field in a three-dimensional semi-infinite composite consists of a Bloch wave, an extinction wave, and an additional term, which is, under certain conditions [stated before Eq. (7) below], spatially localized near the interface. This term can be interpreted as an SPP or as a surface guided wave. Let us assume for simplicity that the composite in question consists of inclusions periodically arranged in the half-space *z* > 0. An incident wave would induce a polarization field **P**(**r**) inside the inclusions (**P** = 0 in vacuum). It is known that the normal modes in an *infinite* composite have the form of the Bloch waves. This applies to every Bloch-periodic function, and in particular to **P**(**r**). However, the Bloch wave is not the normal modes of a semi-infinite composite. Therefore, the polarization induced in such a composite is a superposition of the Bloch wave **P*** _{B}*(

**r**) and an additional surface wave

**P**

*(*

_{S}**r**).

The first of these functions, as any Bloch wave, can be obtained by considering an infinite periodic medium. But the second function is, strictly speaking, not a Bloch wave (it can be viewed as a superposition of infinitely many evanescent Bloch waves [10]), since it does not obey Bloch periodicity in the depth direction. It must be computed by considering the surface explicitly, from the surface-wave equation, which, in the frequency domain, is of the form [21] (here we use notations slightly different from those used in [21]):

**r**∈ Ω

_{tot}, where the latter is the spatial region occupied by all inclusions,

*χ*= (3/4

*π*)(

*ε*− 1)/(

*ε*+ 2) is the coupling parameter,

*ε*is the dielectric permittivity of the inclusions at the working frequency, and

**F**

*(*

_{S}**r**) is a free term, which is completely determined by the transmitted Bloch wave. The function

**F**

*(*

_{S}**r**) was computed explicitly in Ref. [21]. We find it useful to adduce this result here in a self-consistent and succinct, yet somewhat less general from. Consider an incident plane wave

**k**

*=*

_{i}**x̂**

*k*+

_{ix}**ẑ**

*k*is the incident wave vector which satisfies ${k}_{i}^{2}={k}_{ix}^{2}+{k}_{iz}^{2}={k}_{0}^{2}={\left(\omega /c\right)}^{2}$. Thus, referring to a rectangular reference plane

_{iz}*XYZ*, the plane of incidence is

*y*= 0 and the interface is located in the plane

*z*= 0. Then

**p**= (2

*π*/

*h*)(

**x̂**

*n*+

_{x}**ŷ**

*n*) (

_{y}*n*,

_{x}*n*= 0, ±1, ±2,...) are the reciprocal lattice vectors with zero component along the

_{y}*z*-axis, ⊗ denotes tensor product,

*q*is the

_{z}*z*-component of the transmitted Bloch wave vector (the

*x*-component is the same as in the incident wave, so that we can write

**q**=

**x̂**

*k*+

_{iz}**ẑ**

*q*), and

_{z}**P**

^{˜}

*(*

_{B}**k**) is the Fourier transform of the transmitted Bloch wave of polarization

**P**

*(*

_{B}**r**) taken over an elementary cell. The transform is defined as follows. Consider an elementary cubic cell

*C*centered at the point

_{n}**r**

*. Then for*

_{n}**r**∈

*C*, we can write

_{n}**P**

*(*

_{B}**r**) = exp(

*i*

**q**·

**r**

*)*

_{n}**P**

_{cell}(

**r**

*+*

_{n}**R**), where the Cartesian components of

**R**satisfy −

*h*/2 ≤

*R*≤

_{α}*h*/2. The Fourier transform is defined as

Equations (1)–(6) completely characterize the surface wave from the mathematical standpoint. However, direct analysis of these equations is complicated and requires, at the very least, finding the Bloch wave vector **q** and the corresponding function **P**_{cell}(**R**) in an absorbing medium. This is conceptually different from the similar problem encountered in conventional photonic crystals composed of transparent and nondispersive dielectrics; in the latter case, the problem can be reduced to a real symmetric eigenproblem whilst the problem encountered here is not even Hermitian. In addition to computing the Bloch wave, one also faces the problem of solving the integral equation Eq. (1), which can have resonant and therefore non-perturbative solutions. Thus, instead of computing the surface wave directly, we have adopted in this paper an approach based on computing the total field and polarization inside the composite, which includes both the direct and reflected Bloch waves (in the case of a finite slab), as well as the surface wave contributions.

Before proceeding, it is useful to point out one aspect of the surface wave that is easily amenable to analysis. If
${\left(\widehat{\mathbf{x}}{k}_{ix}+\mathbf{p}\right)}^{2}>{k}_{0}^{2}\forall \mathbf{p}\ne 0$, then all the quantities *Q*** _{p}** in (3) have nonzero imaginary parts. Therefore, the function

**F**

*(*

_{S}**r**) decays exponentially away from the interface. The same is true in this case for the surface wave of polarization,

**P**

*(*

_{S}**r**). In the homogenization limit, the exponential decay is fast. Indeed, in the limit

*h*→ 0, we have (for

**p**≠ 0):

*Q*

**→**

_{p}*i*|

**p**|,

**k**→

_{p}**p**+

*i*

**ẑ**|

**p**|, exp[

*i*(

*Q*

**−**

_{p}*q*)

_{z}*h*] − 1 → −1|

**g**|

*h*. With these limits taken into account, the surface wave takes the following form:

**E**

*(*

_{S}**r**) decays exponentially on the scale of

*h*and is, therefore, evanescent, as any SPP.

## 3. Effective medium theory

In this paper, we investigate the influence of the surface wave on the effective medium parameters. To this end, we use the non-asymptotic homogenization theory of Refs. [15–17]. This theory uses accurate approximations of the exact electromagnetic fields **b**, **h**, **e** and **d** in the composite to define the coarse-grained fields **B**, **H**, **E**, **D**. It should be emphasized that by **b**, **h**, **e**, and **d** are not “truly microscopic”, e.g., atomic-scale fields, but rather rapidly varying fields on the scale of tens of nanometers. These fields still obey macroscopic Maxwell’s equations in the composite with spatially-varying material parameters. In contrast, the “macroscopic” or, as we refer to them, coarse-grained fields **B**, **H**, **E**, **D** are obtained by appropriate interpolation procedures and experience spatial variations on larger spatial scales. Note that in intrinsically-nonmagnetic composites, **b** = **h** identically.

The key premise of the homogenization method of Refs. [15–17] is that the coarse-grained fields must satisfy Maxwell’s equations everywhere in space including the sample boundary. Consequently, **E** and **H** are sought as *curl-conforming* interpolants of the rapidly varying fields **e**, **b** – i.e., as the interpolants which preserve the tangential continuity across all interfaces. At the same time, **B**, **D** are sought as the *div-conforming* interpolants which preserve the normal continuity across all interfaces. The procedure is closely related to the theory of discrete Hodge operators [22, 23].

Once the coarse-grained fields have been computed, a linear map *ℒ* is sought between the field pairs (**E**, **H**) and (**D**, **B**). In a suitable “canonical” basis, the operator *ℒ* becomes a generalized material tensor with a leading 6 × 6 block relating (**D**, **B**) to (**E**, **H**). Note that *ℒ* can also contain a “nonlocal” block that relates the coarse-grained fields to field variations over a cell [17].

The fields **e**, **d**, **h** and **b** in the composite can be approximated by a superposition of suitable basis functions [15, 17].

The algorithmic steps for obtaining the generalized material tensor are as follows [15–17]:

- Choose a set of
*N*approximating modes. - Choose a set of
*M*and^{EH}*M*degrees of freedom (d.o.f.s) for the (^{DB}**E**,**H**) and (**D**,**B**) pairs, respectively. The d.o.f. will in general include the mean values of the tangential components of**E**,**H**and of the normal components of**D**,**B**; in addition, the mean values of some derivatives of**E**,**H**may be included. By increasing the number of d.o.f., one trades higher accuracy for a greater level of nonlocality in the characterization of the material. Typically for 3D problems,*M*= 6 (three mean values for each of the two fields) but^{DB}*M*≥ 6. Note that nonlocal d.o.f.s may be included in addition to the mean values.^{EH} - For each mode
*m*= 1, 2,...*M*, compute its respective d.o.f. (the mean boundary values of the tangential components of**E**,**H**for this mode, etc.) Assemble the d.o.f. for the**E**,**H**fields into the*m*th column of matrix*W*and the d.o.f. for the^{EH}**D**,**B**fields into the*m*th column of matrix*W*. Ultimately, matrix^{DB}*W*is of dimension^{EH}*M*×^{EH}*N*and matrix*W*is^{DB}*M*×^{DB}*N*(typically 3 ×*N*).

The procedure employed in this paper follows [15] and involves, as an additional step, volume averaging of auxiliary material tensors defined pointwise, as explained in [15]. The approximations involved in this procedure and the errors incurred are described in detail in [15–17]. In particular, the “in-the-basis error” *γ* that comes from the least squares fit of the material relation (8) is defined as

## 4. Numerical results

We consider a composite consisting of periodically-arranged cubic cells of a period *h*. A spherical gold particle of radius *a* < *h*/2 is located at the center of each cell; the rest of the cell is assumed to be air. The sample under investigation is a plane-parallel slab that contains a finite number *N* of cells in the *z*-direction but is infinite in the *x*- and *y*-directions. The crystallographic axes of the composite coincide with the axes of the laboratory frame *xyz*.

Simulations were performed using the public-domain FDTD package MEEP [24]. Since we consider the case when the slab is illuminated by an infinite-front plane wave that is incident on the slab in the +*z* direction (Fig. 1), the electromagnetic fields inside the composite satisfy Bloch-periodic boundary conditions. Therefore, the computational domain is effectively reduced to *N* elementary cells arranged as shown in Fig. 1(a). In all simulations reported below, the radius of the sphere is *a* = 20 nm, the lattice unit is *h* = 80 nm, so that the filling factor is *f* = 4*πa*^{3}/3*h*^{3} ≈ 0.07. We consider the spectral interval 300 nm ≤ *λ* ≤ 900 nm. The unit of the FDTD cubic grid is Δ = 2.857 nm. In MEEP, boundary conditions at the exterior boundary are enforced by adding perfectly matched layers (PMLs). Additionally, MEEP allows one to apply Bloch-periodic boundary conditions in the *x*- and *y*-directions. The Lorentz-Drude model with six pole expansion terms is used for modeling gold particles. The material parameters are taken from [25].

The basis set used in this paper contained twelve functions (6 incidence directions, ±**x̂**, ±**ŷ**, and ±**ẑ** and two transverse polarizations for each incidence direction). One can, of course, take advantage of the symmetries of the lattice cell and reduce the number of the actual FDTD simulations. For a cubic lattice of spheres, only one basis function needs to be computed numerically for a cell in the bulk and four functions for a cell at the surface; other basis functions can be obtained just by symmetry and rotation. In general, the number of basis functions needed to obtain an accurate homogenization result and to quantify the errors incurred may be greater than 12.

The total field is a superposition of a surface wave and Bloch waves. To picture the surface wave – or equivalently, to remove the Bloch components from the total field – we follow a procedure similar to [26]. More specifically, we considered 13 layers of gold particles with the geometric parameters specified above and found the best fit of the form *E*_{0} exp(*iK _{B}z_{n}*) to the field values

*E*at the cell boundaries

_{y}*z*

_{1},

*z*

_{2}, …

*z*

_{7}of the seven inner layers. Color plots of the total field, its Bloch component and the surface wave are pictured in Fig. 2 for

*λ*= 20

*h*/3. The slab is illuminated with a plane wave propagating along the +

*z*direction. As seen in Fig. 2(c), the surface wave decays rapidly in the direction normal to the surface.

We next consider a composite consisting of *L* = 5 layers and compare the effective parameters of a cell adjacent to the surface (layer *l* = 1; cyan circles in Fig. 3) with those of a cell in the center of the slab (*l* = 3; green diamonds in Fig. 3). An appreciable difference between the two cases is observed, which indicates a noticeable effect of the surface wave in this case. The ”in-the-basis error” *γ* (9) is plotted in Fig. 3(d).

Figure 3 also addresses a somewhat different question. Namely, we assume that the “central” cell is used to compute the effective parameters and wish to estimate the total number of layers for which the parameters obtained in this manner accurately represent the bulk material. In other words, we want to know the minimum value of *L* for which the effective parameters computed using the “central” cell no longer depend on *L*. The effective parameters shown in Fig. 3 for *L* = 5 (green diamonds) and *L* = 9 (red triangles) are also compared to those obtained from Lewin’s theory [27] (solid blue line). The discrepancy between the *L* = 5 and *L* = 9 cases is small, which indicates that, for the parameters considered, *L* = 5 is sufficient to represent bulk samples.

Applicability of the effective parameter description to predicting the transmission (*T*) and reflection (*R*) coefficients of a composite sample is illustrated in Fig. 4 where we compare *T* and *R* of a homogenized slab to the “brute force” FDTD simulations in the composite structure with spherical inclusions present. We used the effective parameters obtained by our theory and also by Lewin’s theory (see plots of the respective effective parameters in Fig. 3). When our homogenization procedures are used, the overall agreement is quite good (this has already been observed in [16]). These results were obtained with bulk parameters only, which indicates that the surface wave does not significantly affect transmission and reflection in the case considered; however, this may not be true in general, e.g., in the case of curved boundaries. Lewin’s theory [27] predicts significantly smaller magnetic effects compared to our homogenization result. More specifically, for small vacuum wave numbers *k*_{0}, Lewin’s prediction for the permeability deviates from unity asymptotically as
${f}_{v}{k}_{0}^{2}{n}_{p}^{2}{a}^{2}/10$, where *f _{v}* and

*n*are the volume fraction and the index of refraction of the particles. This correction tends to be relatively small for the parameters considered, as was already noted [16]. For the absolute values of the transmission and reflection coefficients, Lewin’s results are less accurate than ours (Fig. 4, top panel). However, Lewin’s parameters predict the phases of the TR coefficients quite well even for shorter wavelengths (Fig. 4, bottom panel). We conjecture that Lewin’s theory captures correctly the first nonvanishing correction to the refractive index of the medium but not necessarily to the impedance. Then the above result can be explained by noting that the phases of

_{p}*R*and, especially, of

*T*can be very sensitive to the refractive index but can in some cases tolerate a moderate error in the impedance.

## 5. Summary and discussion

Our homogenization procedure [15–17] is based on the analysis of electromagnetic fields inside the composite and, as such, allows one to distinguish between the corresponding effective parameters in the bulk and near the surface. We anticipate that application of these methods will prove instrumental in developing accurate homogenization models needed to design structures and devices with periodic electromagnetic composites. Although the qualitative difference between surface and bulk parameters has been previously noted [2, 4, 8, 9, 14] and could have been reasonably anticipated, the methods developed in the present paper can be used for rigorous quantitative analysis of the problem.

Surface effects and waves are quite complex and material-dependent, but the paper shows that it is possible to rigorously and quantitatively distinguish between the surface and bulk parameters, or even introduce position-dependent parameters. This approach to homogenization allows one to improve the accuracy of analysis in two different ways. First, under the constraints of *local* approximations of parameters, one can solve Maxwell’s equations numerically using standard methods but with position-dependent parameter tensors. Secondly, extended material tensors that represent nonlocal effects [17] could be used. While the onset of spatial dispersion (nonlocality) is expected to be detrimental for some of the applications that are currently being discussed in the literature (most notably, super-resolution and cloaking), numerical consideration of nonlocality may be worthwhile for other applications of electromagnetic composites, particularly in the situations where full-scale simulations of the microstructure are prohibitively expensive. We note, however, that the apparent “nonlocal” behavior or the constitutive parameters that arises in the theory of homogenization should not be viewed as a direct or complete physical analog of the well-known nonlocal effects in natural materials, such as optical activity in sugar solutions, etc. The reason is that physical nonlocality is truly a microscopic effect that cannot be understood within the macroscopic Maxwell theory. In contrast, fields in metamaterials (with given intrinsic parameters of its constituents) are fully describable by Maxwell’s equations.

## Acknowledgments

This research was supported in part by the Research Grants Council of Hong Kong ( GRF 713011 and GRF 712612), National Science Foundation of China ( NSFC 61271158), HKU 201102160033, and by the University Grants Council of Hong Kong ( AoE/P-04/08). IT and VM acknowledge support from the US National Science Foundation under Grant DMS1216970.

## References and links

**1. **J. D. Baena, L. Jelinek, R. Marques, and M. Silveirinha, “Unified homogenization theory for magnetoinductive and electromagnetic waves in split-ring metamaterials,” Phys. Rev. A **78**, 013842 (2008) [CrossRef] .

**2. **C. R. Simovski, “Material paremeters of metamaterials (a review),” Opt. Spectrosc. **107**, 766–793 (2009) [CrossRef] .

**3. **C. Fietz and G. Shvets, “Current-driven metamaterial homogenization,” Physica B **405**, 2930–2934 (2010) [CrossRef] .

**4. **C. R. Simovski, “On electromagnetic characterization and homogenization of nanostructured metamaterials,” J. Opt. **13**, 103001 (2011) [CrossRef] .

**5. **R. V. Craster, J. Kaplunov, E. Nolde, and S. Guenneau, “High-frequency homogenization for checkerboard structures: defect modes, ultrarefraction, and all-angle negative refraction,” J. Opt. Soc. Am. A **28**, 1032–1040 (2011) [CrossRef] .

**6. **S. Guenneau and F. Zolla, “Homogenization of three-dimensional finite photonic crystals,” Prog. Electromagnetic Res. **27**, 91–127 (2011) [CrossRef] .

**7. **Y. Liu, S. Guenneau, and B. Gralak, “A route to all frequency homogenization of periodic structures,” (2012). ArXiv:1210.6171v2.

**8. **C. Menzel, T. Paul, C. Rockstuhl, T. Pertsch, S. Tretyakov, and F. Lederer, “Validity of effective material parameters for optical fishnet metamaterials,” Phys. Rev. B **81**, 035320 (2010) [CrossRef] .

**9. **C. Menzel, C. Rockstuhl, R. Iliew, F. Lederer, A. Andryieuski, R. Malureanu, and A. V. Lavrinenko, “High symmetry versus optical isotropy of a negative-index metamaterial,” Phys. Rev. B **81**, 195123 (2010) [CrossRef] .

**10. **T. Paul, C. Menzel, W. Smigaj, C. Rockstuhl, P. Lalanne, and F. Lederer, “Reflection and transmission of light at periodic layered metamaterial films,” Phys. Rev. B **84**, 115142 (2011) [CrossRef] .

**11. **C. Rockstuhl, T. Paul, F. Lederer, T. Pertsch, T. Zentgraf, T. P. Meyrath, and H. Giessen, “Transition from thin-film to bulk properties of metamaterials,” Phys. Rev. B **77**, 035126 (2008) [CrossRef] .

**12. **D. R. Smith, S. Schultz, P. Markoš, and C. M. Soukoulis, “Determination of effective permittivity and permeability of metamaterials from reflection and transmission coefficients,” Phys. Rev. B **65**, 195104 (2002) [CrossRef] .

**13. **X. Chen, B. I. Wu, J. A. Kong, and T. M. Grzegorczyk, “Retrieval of the effective constitutive parameters of bianisotropic metamaterials,” Phys. Rev. E **71**, 046610 (2005) [CrossRef] .

**14. **C. Menzel, C. Rockstuhl, T. Paul, F. Lederer, and T. Pertsch, “Retrieving effective parameters for metamaterials at oblique incidence,” Phys. Rev. B **77**, 195328 (2008) [CrossRef] .

**15. **I. Tsukerman, “Effective parameters of metamaterials: a rigorous homogenization theory via Whitney interpolation,” J. Opt. Soc. Am. B **28**, 577–586 (2011) [CrossRef] .

**16. **A. Pors, I. Tsukerman, and S. I. Bozhevolnyi, “Effective constitutive parameters of plasmonic metamaterials: Homogenization by dual field interpolation,” Phys. Rev. E **84**, 016609 (2011) [CrossRef] .

**17. **I. Tsukerman, “Nonlocal homogenization of metamaterials by dual interpolation of fields,” J. Opt. Soc. Am. B **28**, 2956–2965 (2011) [CrossRef] .

**18. **V. A. Markel and J. C. Schotland, “On the sign of refraction in anisotropic non-magnetic media,” J. Opt. **12**, 015104 (2010) [CrossRef] .

**19. **H. Raether, *Surface Plasmons on Smooth and Rough Surfaces and on Gratings* (Springer-Verlag, New York, 1988).

**20. **P. A. Belov and C. R. Simovski, “Boundary conditions for interfaces of electromagnetic crystals and the generalized Ewald-Oseen extinction principle,” Phys. Rev. B **73**, 045102 (2006) [CrossRef] .

**21. **V. A. Markel and J. C. Schotland, “Homogenization of Maxwell’s equations in periodic composites,” Phys. Rev. E **85**, 066603 (2012) [CrossRef] .

**22. **R. Hiptmair, “Discrete Hodge operators,” Num. Math. **90**, 265 (2001) [CrossRef] .

**23. **R. Hiptmair and I. Tsukerman, “Non-asymptotic homogenization of electromagnetic metamaterials via discrete Hodge operators with Trefftz calibration,” to appear in *COMPUMAG – Proceedings of 19th International Conference on the Computation of Electromagnetic Fields* (Budapest, 2013).

**24. **A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. Johnson, “MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method,” Comp. Phys. Comm. **181**, 687 (2010) [CrossRef] .

**25. **A. D. Rakic, A. B. Djurisic, J. M. Elazar, and M. L. Majewski, “Optical properties of metallic films for vertical-cavity optoelectronic devices,” Appl. Opt. **37**, 5271 (1998) [CrossRef] .

**26. **S. Ha, A. A. Sukhorukov, D. K. B., L. C. Botten, C. M. de Sterke, and Y. S. Kivshar, “Bloch-mode extraction from near-field data in periodic waveguides,” Opt. Lett. **34**, 3776 (2009) [CrossRef] [PubMed] .

**27. **L. Lewin, “The electrical constants of a material loaded with spherical particles,” Proc. Inst. Elec. Eng. **94**, 65–68 (1947).