## Abstract

We present a very efficient recursive method to calculate the effective optical response of metamaterials made up of arbitrarily shaped inclusions arranged in periodic 3D arrays. We apply it to dielectric particles embedded in a metal matrix with a lattice constant much smaller than the wavelength of the incident field, so that we may neglect retardation and factor the geometrical properties from the properties of the materials. If the conducting phase is continuous the low frequency behavior is metallic, and if the conducting paths are thin, the high frequency behavior is dielectric. Thus, extraordinary-transparency bands may develop at intermediate frequencies, whose properties may be tuned by geometrical manipulation.

© 2010 Optical Society of America

## 1. Introduction

Metallic films with sub-wavelength nanometric holes may display an extraordinarily large transmittance at near infrared frequencies for which the metal is opaque and light waves are not expected to propagate within the holes [1, 2]. The anomalous transmission and other extraordinary optical properties of nano-structured metallic films open the possibility of tailored design for many applications that include hyperlens far-field-subdiffraction imaging [3–5], cloaking [6], optical antennas [7, 8], and circular polarizers [9]. Thus, understanding the electromagnetic properties of this sometimes called plasmonic metamaterials has become important.

Many different approaches have been proposed to describe the extraordinary optical properties of some structured systems [10–26]. Some resonances have been identified with surface-plasmon-polaritons (SPP’s) that may be excited by light after being scattered by the system.

In a recent work [27] the frequency-dependent complex macroscopic dielectric-response tensor ${\varepsilon}_{ij}^{M}(\omega )$ of 2D-periodic lattices of cylindrical inclusions with arbitrarily shaped cross sections embedded within metallic hosts were obtained in the local limit and enhanced transparency was obtained without invoking explicitly a SPP mechanism. In this paper we develop Haydock’s recursive scheme [28] to obtain the macroscopic dielectric response of 2D and 3D periodic metamaterials in the long wavelength limit. Our method yields a speed improvement of several orders of magnitude over that of Ref. [27], which allows previously untractable calculations for 3D structures with arbitrary geometry, such as interpenetrated inclusions made out of dispersive and dissipative components. We show that the geometry of the inclusions and of the lattice might lead to very anisotropic optical behavior and to a very generic enhanced transmittance for metal-dielectric metamaterials whenever there are only poor conducting paths across the whole sample.

## 2. Theory

We consider a periodic lattice of arbitrarily shaped nanometric inclusions (*b*) embedded within a homogeneous material (*a*). We assume that each region *α* = *a, b* is large enough to have a well defined macroscopic response *ɛ _{α}* which we assume local and isotropic, but much smaller than the free wavelength

*λ*

_{0}= 2

*πc*/

*ω*with

*c*the speed of light in vacuum and

*ω*the frequency. The microscopic response is described by

*ɛ*≡

_{ab}*ɛ*–

_{a}*ɛ*and

_{b}*B*(

**) = 0 or 1 is the characteristic function for the**

*r**b*regions, which we assume periodic,

*B*(

**+**

*r***) =**

*R**B*(

**), with {**

*r***} the Bravais lattice of the metamaterial. The constitutive equation**

*R***(**

*D***) =**

*r**ɛ*(

**)**

*r***(**

*E***) may be written in reciprocal space as**

*r***(**

*D***) is the displacement field and**

*r***(**

*E***) the electric field,**

*r***(**

*D*_{G}**) and**

*q***(**

*E*_{G}**) the corresponding Fourier coefficients with wavevector**

*q***+**

*q***,**

*G***is the conserved Bloch’s vector and {**

*q***} the reciprocal lattice. Here,**

*G***is the**

*ɛ*_{GG′}**–**

*G***Fourier coefficient of**

*G′**ɛ*(

**). We consider now a longitudinal external field**

*r*

*E*^{ex}(

**) = −∇**

*r**ϕ*

^{ex}(

**) and we neglect retardation within the small unit cell, so we may assume that the total electric field is**

*r**longitudinal*

**+**

*q***)/|**

*G***+**

*q***| by**

*G***. We remark that our final results apply as well to**

*Ĝ**transverse*fields. As ∇ ×

*E*^{ex}= ∇ ×

*D**=*

^{L}**0**and ∇ ·

*E*^{ex}= ∇ ·

*D**= 4*

^{L}*πρ*

^{ex}, we may identify

*E*^{ex}with

*D**which we chose as a plane wave with wavevector*

^{L}**without small lengthscale fluctuations, ${\mathit{\text{D}}}_{\mathit{\text{G}}\ne 0}^{L}(\mathit{\text{q}})=0$. Substituting Eq. (3) into the longitudinal projection of Eq. (2) allows us to solve for**

*q***00**component. The macroscopic longitudinal field ${\mathit{\text{E}}}_{M}^{L}$ is obtained from

*E**by eliminating its spatial fluctuations, i.e., ${\mathit{\text{E}}}_{M}^{L}={\mathit{\text{E}}}_{0}^{L}$. Similarly, ${\mathit{\text{D}}}_{M}^{L}={\mathit{\text{D}}}_{0}^{L}$. Thus, from Eq. (4) we identify the longitudinal projection of the macroscopic dielectric response*

^{L}We should emphasize that *ξ* in Eq. (6) depends in general on the direction of ** q**. Calculating
${\mathit{\text{\u025b}}}_{ML}^{-1}=\widehat{\mathit{\text{q}}}\widehat{\mathit{\text{q}}}\cdot {\mathit{\text{\u025b}}}_{M}^{-1}\cdot \widehat{\mathit{\text{q}}}\widehat{\mathit{\text{q}}}$ for several propagation directions

**[32] we may obtain**

*q̂**all*the components of the

*full*inverse long-wavelength dielectric tensor ${\mathit{\text{\u025b}}}_{M}^{-1}$ and from it

*ɛ*_{M}. Properties such as reflectance and transmittance may then be calculated using standard formulae [33, 34].

Equation (6) is the starting point for a very efficient numerical calculation of the macroscopic response of nanostructured metamaterials, as developed in the next section.

## 3. Haydock’s recursive method

To calculate efficiently the macroscopic response, we start by taking the Fourier transform of Eq. (1), *ɛ _{GG′}* =

*ɛ*–

_{a}δ_{GG′}*ɛ*, where

_{ab}B_{GG′}*v*within the unit cell of volume Ω. In particular,

*B*

**=**

_{00}*v*/Ω ≡

*f*is the filling fraction of the inclusions. From Eq. (5) we obtain

*𝒢*(

_{GG′}*u*) are the matrix elements of a Green’s operator

*𝒢̂*, i.e., the resolvent of the operator (

*u*–

*ℋ̂*), corresponding to a

*Hamiltonian ℋ̂*with elements

*spectral variable u*(

*ω*) ≡ (1 –

*ɛ*(

_{b}*ω*)/

*ɛ*(ω))

_{a}^{−1}plays the role of

*energy*.

From Eqs. (6) and (8) we obtain *ξ* = *𝒢*** _{00}**(

*u*)/

*ɛ*= 〈

_{ab}**0**|

*𝒢̂*(

*u*)|

**0**〉/

*ɛ*, where |

_{ab}**0**〉 represents the state corresponding to a plane wave with wave vector

**, and we denote in general the state corresponding to a plane wave with wave vector**

*q***+**

*q***as |**

*G***〉. This allows the use of Haydock’s recursive scheme [28] to obtain the projected Green’s function and thus the macroscopic response. We set |−1〉 = 0, |0〉 = |**

*G***0**〉,

*b*

_{0}= 0 and recursively define the states |

*n*〉 through

*a*and

_{n}*b*can be recursively obtained by demanding that the states |

_{n}*n*〉 be normalized and orthogonal to each other, yielding

*a*

_{n}_{−1}= 〈

*n*– 1|

*ñ*〉 = 〈

*n*– 1|

*ℋ̂*|

*n*– 1〉 and ${b}_{n}^{2}=\u3008\tilde{n}|\tilde{n}\u3009-{a}_{n-1}^{2}-{b}_{n-1}^{2}$.

From Eq. (9) we notice that the action of *ℋ̂* includes first of all a multiplication with a vector (in cartesian space) operator with elements *Ĝ**δ _{GG′}*, which is trivially done in

*reciprocal*space. Afterwards we have to multiply with a scalar operator with elements

*B*. However, from Eq. (7), this product is actually a convolution. Thus, a fast Fourier transform allows us to perform this operator in

_{GG′}*real*space simply as a multiplication with the characteristic function

*B*(

**). Finally, we transform back into reciprocal space and apply an inner product with**

*r*

*Ĝ**δ*. This way, by going back and forth between real and reciprocal space we can calculate the action of

_{GG′}*ℋ̂*performing only trivial multiplications with diagonal matrices. We implemented the calculation of Haydock’s coefficients using the

*Perl Data Language*[35, 36].

In the basis {|*n*〉} the Hamiltonian is represented by a symmetric tridiagonal matrix with main diagonal elements {*a _{n}*} and sub/sup-diagonal elements {

*b*}. Then, we may write

_{n}*𝒢̂*

^{–1}=

*u*–

*ℋ̂*recursively in blocks as

*A*= (

_{n}*u*–

*a*) and

_{n}*ℬ*= (−

_{n}*b*, 0, 0, ⋯), and where

_{n}*𝒢̂*≡

*𝒢̂*

_{0}. Here we used calligraphic letters to denote any matrix except 1 × 1 matrices which are equivalent to scalars.

Now we write *𝒢 _{n}* in blocks as

*𝒢*

**(**

_{00}*u*) =

*R*

_{0}given recursively through the continued fraction

Notice that Haydock’s coefficients depend only on the geometry through
${B}_{\mathit{\text{GG}}\prime}^{L}$. The dependence on composition and frequency is completely encoded in the complex valued spectral variable *u*. Thus, for a given geometry we may explore manifold compositions and frequencies without recalculating Haydock’s coefficients. Furthermore, our method is equally suited to the calculation of the properties of systems made up of non-dispersive, transparent materials as well as dispersive, opaque and dissipative media; we simply have to provide the appropriate (complex) value for *u*.

## 4. Results

To test our calculation scheme, in Fig. 1(a) we show the macroscopic response *ɛ _{M}* calculated for a simple cubic lattice of dielectric spheres, with dielectric function

*ɛ*= 4, embedded within a dielectrically inert host

_{b}*ɛ*= 1 for various filling fractions

_{a}*f*. We compare the results of our formulation based on Haydock’s recursive scheme (H) with the standard Maxwell-Garnett (MG) and Bruggeman (B) effective medium theories for 3D. As expected, the three theories coincide at low filling fractions, but not surprisingly, they differ as

*f*increases, with our results being between MG which lies below and B which lie above. Bruggeman theory treats both

*a*and

*b*media on an equivalent basis, while medium

*a*plays the role of the host and medium

*b*that of inclusions in both Maxwell-Garnett theory and our’s. Thus the MG results are closer to ours than B. Similarly, in Fig. 1(b) we show the response

*ɛ*normal to the optical axis of a system composed of a 2D square array of square prisms whose diagonals are oriented along the conventional unit cell sides, i.e., the prisms are rotated 45° with respect to the unit cell. As in the 3D case, we took

_{M}*ɛ*= 4 and

_{b}*ɛ*= 1. We compare our theory with the 2D MG and B formulae. As in the 3D example, the three formulations coincide in the dilute

_{a}*f*→ 0 limit. For larger

*f*, MG yields the lowest

*ɛ*and

_{M}*B*the largest. At the percolation limit

*f*= 0.5, at which the edges of nearby prisms touch each other, the geometry does not distinguish between inclusions and host. At this limit, our results coincide with the symmetric Bruggeman 2D formula. Thus, in this particular 2D case, B lies closer to our results than MG. Furthermore, due to the symmetry of this system at

*f*= 0.5, Keller’s theorem [38–40] implies that ${\varepsilon}_{M}=\sqrt{{\varepsilon}_{a}{\varepsilon}_{b}}$, which is exactly fulfiled by both B and H theories.

In a recent paper [27] an alternative, essentially exact but much more numerically expensive method to calculate the macroscopic response of metamaterials was developed and it was applied to manifold 2D systems. The results were satisfactorily tested against previous calculations and validated through the fulfilment of Keller’s reciprocal theorem [38–40]. We have repeated *all* the calculations in that paper and have verified that our results agree in the long wavelength limit and are about four orders of magnitude faster. The speed enhancement is much larger in 3D. Further 2D calculations with the present formulation were presented in Ref. [37].

In Fig. 2(a) we show the normal-incidence reflectance of an isotropic semi-infinite metamaterial consisting of a simple cubic lattice of small cubical voids embedded with filling fraction *f* = 0.6 within a Drude conductor with plasma frequency *ω _{p}* and a low dissipation parameter Γ = 0.01

*ω*. The low frequency reflectance is almost unity but attains a deep minimum at

_{p}*ω*≈ 0.51

*ω*and becomes very small for 0.74 <

_{p}*ω*/

*ω*< 0.91, it attains a peak at 0.97

_{p}*ω*and finally decreases at high frequencies. This behavior may be understood by looking at Fig. 2(b), which shows the corresponding macroscopic dielectric function. We remark that at low frequencies

_{p}*ɛ*is negative and large corresponding to a metallic behavior, as the interstitial conducting phase percolates. Nevertheless, at larger frequencies there is a series of resonances originated from the coupled plasmon modes of the cubic inclusions [41]. Thus, the macroscopic response becomes dielectric like and approaches at some intermediate frequencies the vacuum value

_{M}*ɛ*= 1. In the figure we have indicated the frequencies

_{V}*ω*≈ 0.80

*ω*and 0.86

_{p}*ω*for which Re[

_{p}*ɛ*] = 1 and the metamaterial would become perfectly transparent were it not for its small Im[

_{M}*ɛ*]. Similarly, at

_{M}*ω*≈ 0.52

*ω*there is a resonance below which Re[

_{p}*ɛ*] approaches 1, corresponding to a sharp minimum in the reflectance followed by a sharp maximum due to resonant absorption. Above the last resonance, at 0.94

_{M}*ω*, there is a narrow frequency band where Re[

_{p}*ɛ*] < 0, corresponding to a maximum in

_{M}*R*. Finally, asymptotically [

*ɛ*] → 1 and

_{M}*R*→ 0.

The figure also shows the transmittance *T* of a *d* = 200 nm wide film made of the same metamaterial. We notice that below *ω _{p}* it is almost null except for a band where it takes appreciable values and which corresponds to the frequencies where the reflectance

*R*of the corresponding semi-infinite system becomes small. The interpretation of the optical properties of a film is encumbered by the interference between the multiply reflected waves and their decay as they cross the film. Thus,

*T*shows a very rich structure originated from resonant oscillations in the imaginary part of $k=(\omega /c)\sqrt{{\varepsilon}_{M}}$ below

*ω*and Fabry-Perot resonances above

_{p}*ω*. Increasing dissipation diminishes both types of oscillations, as illustrated in the figure for Γ = 0.1

_{p}*ω*without eliminating the transmission band. Finally, as a reference we show the transmittance of low and high dissipation homogeneous Drude films of effective width

_{p}*d*

_{eff}= (1 –

*f*)

*d*= 80 nm, so that it contains the same amount of metal as the metamaterial. Notice that for both cases the transmittance of the metamaterial can be enhanced by about two orders of magnitude above that of the effective film.

In Fig. 3(a) we show the transmittance *T* of a film made of a simple cubic lattice of spherical dielectric inclusions with response *ɛ _{b}* = 4 within an Au host [42]. We chose the radius as

*r*= 0.6

*a*, with

*a*the lattice parameter, so the spheres overlap their neighbors. We have normalized the results to the transmittance

*T*

_{eff}of an effective homogeneous Au film of width

*d*

_{eff}, in order to emphasize the transmittance enhancement due to the metamaterial geometry. Several enhancement peaks between one and two orders of magnitude are visible in the transmittance spectrum, corresponding to the excitation of coupled multipolar plasmon resonances within the region where the metal is opaque.

In Fig. 3(b) we show the normalized transmittance *T _{α}*/

*T*

_{eff}(

*α*=

*x, y*) for plane polarized light normally incident on a 200 nm film lying on the

*xy*plane made of a simple orthorhombic lattice of

*z*-oriented dielectric cylinders with radius

*r*= 0.53

*a*, height

_{x}*h*= 0.9

*a*dielectric function

_{z}*ɛ*= 2 and lattice parameters

_{b}*a*= 1.15

_{y}*a*and

_{x}*a*=

_{z}*a*within an Au host. There is a huge anisotropy, with a peak enhancement of

_{x}*T*almost two orders of magnitude larger than that of

_{y}*T*. For this geometry there is an overlap between neighbor cylinders along

_{x}*x*, so the system is a better low frequency conductor along

*y*.

We should discuss some limitations of our theory. Our assumption of periodicity is not too restrictive, as the efficiency of our scheme permits calculations for complex unit cells with many inclusions. Our use of the long-wavelength approximation does impose a bound on the size *a* of the unit cell. Retardation corrections are expected to be of order (*a*/*λ*_{0})2 and we have verified through comparisons to exact calculations in 2D [27] that for *a* ≤ *λ*_{0}/5 they may be considered negligible for many purposes. The interesting regions in Figs. 3(a) and 3(b) lie below ∼ 2 eV, where our macroscopic response would be accurate up to *a* ∼ 100 nm, a scale easily attainable with current fabrication techniques. Nevertheless, our theory does neglect the presence of transition layers of width *d _{t}* ∼

*a*close to the metamaterial boundaries which are not well described by a local, position independent macroscopic response

*ɛ*, as inclusions closer to the surface than their size would be sliced by the surface, and the effect of the surface on the microscopic fluctuations of the field extendes a distance ∼

_{M}*a*. Similar transition regions are present at ordinary crystalline and spatially dispersive materials and are known to produce corrections to the optical properties of semi-infinite media of the order of

*d*/

_{t}*λ*

_{0}[43]. For a free-standing subwavelength thin film of width

*d*, the corresponding corrections can become of the order of

*d*/

_{t}*d*. In our case, both corrections would be of only a few percent if the inclusions were of the order of a few nanometers.

In summary, we developed a formalism that allows very efficient calculations of the macroscopic response of 3D nano-structured periodic metamaterials. We applied it to films made of various lattices of dielectric inclusions with assorted shapes within opaque metallic hosts and found an extraordinary enhancement of the transmittance as a generic property whenever the metamaterial is conducting at low frequencies and has dielectric like resonances at larger frequencies.

## Acknowledgments

This work was supported DGAPA-UNAM (IN120909 (WLM)), CONACyT (48915-F (BMS) and by ANPCyT-UNNE (204 y 190-PICTO-UNNE-2007 (GPO)).

## References and links

**1. **L. Novotny and B. Hecht, *Principles of Nano-Optics* (Cambridge University Press, New York, 2006).

**2. **J. Weiner, “The physics of light transmission through subwavelength apertures and aperture arrays,” Rep. Prog. Phys. **72**, 0644011 (2009).

**3. **L. Chen and G. P. Wang, “Pyramid-shaped hyperlenses for three-dimensional subdiffraction optical imaging,” Opt. Express **17**, 3903 (2009).

**4. **Z. Liu, H. Lee, Y. Xiong, C. Sun, and X. Zhang, “Far-field optical hyperlens magnifying sub-diffraction-limited objects,” Science **315**, 1686 (2007).

**5. **A. Salandrino and N. Engheta, “Far-field subdiffraction optical microscopy using metamaterial crystal: theory and simulations,” Phys. Rev. B **74**, 075103 (2006).

**6. **W. Cai, U. K. Chettiar, A. V. Kildishev, and V. M. Shalaev, “Nonmagnetic cloak with minimized scattering,” Appl. Phys. Lett. **91**, 1111051 (2007).

**7. **L. Novotny, “Effective wavelength scaling for optical antennas,” Phys. Rev. Lett. **98**, 2668021 (2007).

**8. **W. Dickson, G. Wurtz, P. Evans, D. O’Connor, R. Atkinson, R. Pollard, and A. Zayats, “Dielectric-loaded plasmonic nanoantenna arrays: A metamaterial with tunable optical properties,” Phys. Rev. B **76**, 115411 (2007).

**9. **J. K. Gansel, M. Thiel, M. S. Rill, M. Decker, K. Bade, V. Saile, G. von Freymann, S. Linden, and M. Wegener, “Gold helix photonic metamaterial as broadband circular polarizer,” *Science*325, 1513 (2009).

**10. **B. Hou, H. Wen, Y. Leng, and W. Wen, “Enhanced transmission of electromagnetic waves through metamaterials,” Appl. Phys. A **87**, 217 (2007).

**11. **L. Zhou, W. Wen, C. T. Chan, and P. Sheng, “Electromagnetic-wave tunneling through negative-permittivity media with high magnetic fields,” Phys. Rev. Lett. **94**, 2439051 (2005).

**12. **J. B. Pendry, L. Martín-Moreno, and F. J. Garcia-Vidal, “Mimicking surface plasmons with structured surfaces,” Science **305**, 847 (2004).

**13. **P. Sheng, R. Stepleman, and P. Sanda, “Exact eigenfunctions for square-wave gratings: application to difraction and surface-plasmon calculations,” Phys. Rev. B **26**, 2907 (1982).

**14. **H. Lochbihler and R. Depine, “Highly conducting wire gratings in the resonance region,” Appl. Opt. **32**, 3459 (1993).

**15. **H. Lochbihler, “Surface polaritons on gold-wire gratings,” Phys. Rev. B **50**, 4795 (1994).

**16. **H. Ghaemi, T. Thio, D. Grupp, T. Ebbesen, and H. Lezec, “Surface plasmon enhanced optical transmission through subwavelength holes,” Phys. Rev. B **58**, 6779 (1998).

**17. **L. Martín-Moreno, F. J. García-Vidal, H. J. Lezec, K. M. Pellerin, T. Thio, J. B. Pendry, and T. W. Ebbesen, “Theory of extraordinary optical transmission through subwavelength hole arrays,” Phys. Rev. Lett. **86**, 1114 (2001).

**18. **D. Grupp, H. Lezec, T. Ebbesen, K. Pellerin, and T. Thio, “Surface plasmon enhanced optical transmission through subwavelength holes,” Appl. Phys. Lett **77**, 1569 (2000).

**19. **S. Darmanyan and A. Zayats, “Light tunneling via resonant surface plasmon polariton states and the enhanced transmission of periodically nanostructured metal films: an analytical study,” Phys. Rev. B **67**, 035424 (2003).

**20. **J. Porto, F. García-Vidal, and J. Pendry, “Transmission resonances on metallic gratings with narrow slits,” Phys. Rev. Lett. **83**, 2845 (1999).

**21. **Q.-H. Park, J. H. Kang, J. W. Lee, and D. S. Kim, “Effective medium description of plasmonic metamaterials,” Opt. Express **15**, 6994 (2007).

**22. **A. Agrawal, Z. Vardeny, and A. Nahata, “Engineering the dielectric function of plasmonic lattices,” Opt. Express **16**, 9601 (2008).

**23. **Q. Cao and P. Lalanne, “Negative role of surface plasmon in the transmission of metallic gratings with very narrow slits,” Phys. Rev. Lett. **88**, 0574031 (2002).

**24. **M. Treacy, “Dynamical diffraction in metallic optical gratings,” Appl. Phys. Lett. **75**, 606 (1999).

**25. **M. Treacy, “Dynamical diffraction explanation of the anomalous transmission of light through metallic gratings,” Phys. Rev. B **66**, 1951051 (2002).

**26. **E. Popov, M. Neviére, S. Enoch, and R. Reinisch, “Theory of light transmission through subwavelength periodic hole arrays,” Phys. Rev. B **62**, 16100 (2000).

**27. **G. P. Ortiz, B. E. Martínez-Zérega, B. S. Mendoza, and W. Mochán, “Effective dielectric response of metamaterials,” Phys. Rev. B **79**, 245132 (2009).

**28. **R. Haydock, “The recursive solution of the Schrödinger equation,” Solid State Phys. **35**, 215 (1980).

**29. **S. L. Alder, “Quantum theory of the dielectric constant in real solids,” Phys. Rev. **126**, 413 (1962).

**30. **N. Wiser, “Dielectric constant with local field effects included,” Phys. Rev. **129**, 62 (1963).

**31. **W. L. Mochán and R. G. Barrera, “Electromagnetic response of systems with spatial fluctuations. i. general formalism,” Phys. Rev. B **32**, 4984 (1985).

**32. **That the long-wavelength response *ɛ _{M}* is independent of the direction of

**→ 0 allows us to use the results of longitudinal calculations to solve optical (i.e., transverse) problems.**

*q***33. **R. G. Barrera, A. Reyes-Coronado, and A. García-Valenzuela, “Nonlocal nature of the electrodynamic response of colloidal systems,” Phys. Rev. B **75**, 184202 (2007).

**34. **M. Born and E. Wolf, *Principles of Optics* (Cambridge University Press, 1999), 7th ed.

**35. **K. Glazebrook, J. Brinchmann, J. Cerney, C. DeForest, D. Hunt, T. Jenness, T. Luka, R. Schwebel, and C. Soeller, “The perl data language v.2.4.4,” Available from http://pdl.perl.org.

**36. **K. Glazebrook and F. Economou, “Pdl: The perl data language,” Dr. Dobb’s Journal (1997). http://www.ddj.com/184410442.

**37. **E. Cortés, W. L. Mochán, B. S. Mendoza, and G. P. Ortiz, “Optical properties of nano-structured metamaterials,” Phys. Status Solidi B **247**, 2102 (2010).

**38. **J. B. Keller, “Conductivity of a medium containing a dense array of perfectly conducting spheres or cylinders or nonconducting cylinders,” J. Appl. Phys. **34**, 991 (1963).

**39. **J. B. Keller, “A theorem on the conductivity of a composite medium,” J. Math. Phys. **5**, 548 (1964).

**40. **J. Nevard and J. B. Keller, “Reciprocal relations for effective conductivities of anisotropic media,” J. Math. Phys. **26**, 2761 (1985).

**41. **R. Fuchs, “Theory of the optical properties of ionic crystal cubes,” Phys. Rev. B **11**, 1732 (1975).

**42. **P. Johnson and R. Christy, “Optical constant of noble metals,” Phys. Rev. B **6**, 4370 (1972).

**43. **W. L. Mochán and R. G. Barrera, “Intrinsic surface-induced optical anisotropies of cubic crystals: local-field effect,” Phys. Rev. Lett. **55**, 1192 (1985).