Semi-phenomenological effective permittivity approach to metallic periodic structures

A detailed review of the theory of effective permittivity for oneand twodimensional periodic structures shows its limited validity for metal-dielectric structures in the visible and near infra-red if the feature dimensions are comparable with the metal skin depth. We propose a phenomenological correction to the static formulae using a realistic assumption for the electric field behavior inside the metal features. This approach allows us to obtain analytical expressions for the effective permittivity in the case when the electric field is not sufficiently homogeneous within the unit cell of the gratings. A comparison with the numerical results of the Fourier modal method demonstrates the validity of the analytical formulae. Additional study is made on the impedance approximation at the outer boundaries of the periodical structure in order to propose analytical formulae for the reflection coefficient that permits better correspondence with the numerical results. The link between the values of effective permittivity and permeability defined as the ratios between the averaged fields, and the metamaterial permittivity and permeability is discussed. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (050.2065) Effective medium theory; (050.6624) Subwavelength structures; (160.3918) Metamaterials; (050.1950) Diffraction gratings. References and links 1. J. C. Maxwell-Garnett, “Colors in metal glasses and in metallic films,” Philos. Trans. R. Soc. London Ser. A 203(359–371), 385– 420 (1904). 2. I. Barzilay, M. L. Myers, L. B. Cooper, and G. N. Graser, “Mechanical and chemical retention of laboratory cured composite to metal surfaces,” J. Prosthet. Dent. 59(2), 131–137 (1988). 3. P. Vukusic, J. R. Sambles, C. R. Lawrence, and R. J. Wootton, “Quantified interference and diffraction in single Morpho butterfly scales,” Proceedings: Biological Sciences, The Royal Society of London 266, 1403–1411 (1999). 4. H. Tada, S. Mann, I. Miaoulis, and P. Wong, “Effects of a butterfly scale microstructure on the iridescent color observed at different angles,” Opt. Express 5(4), 87–92 (1999). 5. S. Tian and H. Brill Robert, Ancient Glass Research along the Silk Road (World Scientific, 2009). 6. G. W. Milton, The Theory of Composites (Cambridge University, 2002). 7. G. W. Milton and K. Golden, “Representations for the Conductivity Functions of Multicomponent Composites,” Commun. Pure Appl. Math. 43(5), 647–671 (1990). 8. R. Clausius, Die Mechanische Behandlung der Electricität (Vieweg + Teubner Verlag, 1879). 9. O. F. Mossotti, “Discussione analitica sull’influenza che l’azione di un mezzo dielettrico ha sulla distribuzione dell’elettricità alla superficie di più corpi elettrici disseminati in esso,” Memorie di Mathematica e di Fisica della Società Italiana della Scienza Residente in Modena 24, 49–74 (1850). 10. D. Yaghjian, “Electric dyadic Green’s functions in the source region,” Proc. IEEE 68(2), 248–263 (1980). 11. P. Lalanne and D. Lemercier-Lalanne, “On the effective medium theory of subwavelength periodic structures,” J. Mod. Opt. 43(10), 2063–2085 (1996). 12. N. A. Nicorovici and R. C. McPhedran, “Transport properties of arrays of elliptical cylinders,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 54(2), 1945–1957 (1996). 13. F. L. Galeener, “Submicroscopic-Void Resonance: The Effect of Internal Roughness on Optical Absorption, Phys. Rev. Lett. 27, 421-423 (1971), “Erratum,” ibid, p.769 (1971). 14. S. M. Rytov, “Electromagnetic properties of a finely stratified medium,” Sov. Phys. JETP 2, 466–475 (1956). 15. P. Lalanne, “Effective medium theory applied to photonic crystals composed of cubic or square cylinders,” Appl. Opt. 35(27), 5369–5380 (1996). Vol. 26, No. 10 | 14 May 2018 | OPTICS EXPRESS 12813 #325782 https://doi.org/10.1364/OE.26.012813 Journal © 2018 Received 9 Mar 2018; revised 29 Mar 2018; accepted 5 Apr 2018; published 3 May 2018 16. L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A 14(10), 2758–2767 (1997). 17. G. Gao, C. Torres-Verdin, and T. M. Habashy, “Analytical Techniques to evaluate the integrals of 3D and 2D spatial Dyadic Green’s functions,” PIERS 52, 47–80 (2005). 18. D. E. Gray, ed., American Institute of Physics Handbook (McGraw-Hill, 1957), sec.6. 19. S. R. Coriell and J. L. Jackson, “Bounds on transport coefficients of two-phase materials,” J. Appl. Phys. 39(10), 4733–4736 (1968). 20. T. Weiss, G. Granet, N. A. Gippius, S. G. Tikhodeev, and H. Giessen, “Matched coordinates and adaptive spatial resolution in the Fourier modal method,” Opt. Express 17(10), 8051–8061 (2009). 21. L. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, and J. R. Andrewartha, “The dielectric lamellar diffraction grating,” Opt. Acta (Lond.) 28(3), 413–428 (1981). 22. J. Y. Suratteau, M. Cadilhac, and R. Petit, “Sur la détermination numérique des efficacités de certains réseaux diélectriques profonds,” J. Opt. (Paris) 14, 273–288 (1983). 23. K. B. Dossou, C. G. Poulton, and L. C. Botten, “Effective impedance modeling of metamaterial structures,” J. Opt. Soc. Am. A 33(3), 361–372 (2016). 24. 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(19), 195104 (2002). 25. J. Yang, C. Sauvan, T. Paul, C. Rockstuhl, F. Lederer, and P. Lalanne, “Retrieving the effective parameters of metamaterials from the single interface scattering problem,” Appl. Phys. Lett. 97(6), 061102 (2010).


Introduction
Effective medium theory plays a major role in explaining and modeling many problems in physics and chemistry.The fundamental work of Maxwell-Garnet [1] refers to works of Stokes, Lorenz, and Lord Rayleigh in the 19 th century that describe the scattering by spheres, small compared to the wavelength of light.During the 20 th century this approach was developed at first in detail for conductivity and transport phenomena for inclusions of different geometries, later extended to electromagnetic and optical properties, as explained in detail in Sec.2.Since more than 20 years a term 'metamaterial' has emerged to describe structures with artificial properties due to the subwavelength structuring and mixing.Curiously, the first use of the term, as far as we can say, was in the dentistry in 1988 [2], quite soon later picked up by opticians.Another remark concerning metamaterials is that nature and humanity have produced and used metamaterials a long time ago.Nanostructuring of some butterfly wings leads to their pigmentless coloring [3,4], while metal inclusions in glass vases or colored window glasses have been used to obtain coloring since 500-400 BC [5].
Effective permittivity (or/and permeability) theories find various applications for qualitative and possibly quantitative analysis of nanotubes (for example, in night vision masks), color pixel filters for visible and IR cameras, photovoltaic and thermic microcells, plasmonic effects, and resonant filters.The common feature of these structures is their periodicity in one (1D) or two (2D) dimensions, whereas photonic crystals, studied extensively in the 90s and beginning of this century can have three-dimensional (3D) periodicity.
While the theory is quite well established for periodic structures having such small feature dimensions, which is considered as almost homogenized because the electromagnetic field varies very slowly inside (see Sec. 2), there are many still open problems to solve.The first one that is the topic of our study concerns metallic gratings with features (much) smaller than the wavelength in vacuum, but larger than the skin depth of the metal.In that case it is not possible to consider the EM field as almost constant inside.
Another problem that lies out of the scope of this paper deals with aperiodic (random) distribution of inclusion and has been the topic of extensive interest during the 20 th century.An interested reader can find a review by Milton [6].In some cases it is even necessary to average not only the values of the field, but also of the intensity (modulus square) of the electric field (ch.16 of [6], and [7]).
In the paper we use the definition of the effective permittivity as the ratio between the average values of the electric displacement and electric field, which differs from the definition used in the recent metamaterial papers, but there is an explicit relation between them, as discussed in Sec. 5.
In the following we treat several different cases with one-or two-dimensional periodicity.The first class of examples considers features with dimensions much smaller than the skin depth, so that the classical effective medium theory, which is based on the assumption that the electric field is constant within the metal, is valid.The second type of examples treats the opposite case with feature dimensions larger than the skin depth, for which the effective medium approach does not work, and we propose a correction based on the modal method.The reference numerical results are obtained using the rigorous Fourier modal method, also known as RCWA (rigorous coupled-wave approach).
Section 2 presents a short review of the results of the effective medium theory, together with some confirmative examples, and points out the problem for larger features.In Sec. 3 we propose a semi-phenomenological solution to this problem, and demonstrate the existence of another discrepancy with numerical results for the reflectivity of the effective layer due to the inconsistency of the presentation of the surface impedance.An attempt to qualitatively solve this problem is given in Sec. 4. Section 5 gives the link between the effective permittivity and permeability defined in the homogenization approach as the relations between the averaged field, and the metamaterial permittivity and permeability.We discuss why the average effective permeability of non-magnetic materials is equal to the vacuum permeability, whereas its metamaterial counterpart can present magnetic properties, and show how the metamaterial permittivity and permeability can be obtained from the results of Sec. 2 -4.The derivation of the Green tensor approach and its singularities are given in detail in Appendices A and B.

Effective medium permittivity for periodic homogenized structures
We consider a system with 1 or 2-dimensional periodicity made of lamellae or cylindrical objects invariant in vertical direction (Fig. 1).In order to obtain a simple model that enables us to better understand the behavior of these short-pitch structures, we use a homogenization approach based on the Maxwell-Garnet approach [1] and Green tensor singularity analysis.Let us suppose that the unit cell contains two regions 1 and 2, within each of them the permittivity j ε (j = 1, 2) is homogeneous.We assume that the region 2 is included in region 1.For 1D grating this assumption is ambiguous, but the formulas that are obtained are symmetrical in that case.In the case of anisotropic medium 2, 2 ε is a tensor.V is the volume of the entire unit cell, V j (j = 1, 2) are the volumes of each region.The filling ratio is f = V 2 /V.The original approach proposed by Maxwell-Garnet [1] enables us to find the effective permittivity ε eff of the inhomogeneous medium, based on a simple hypothesis that it gives the link between the mean values (arithmetic means) of electric and displacement fields: where a pair of angular brackets stands for the arithmetic mean value and ε 0 is the vacuum permittivity.Similar definition applies for effective permeability, which for non-magnetic media is simply equal to the vacuum permeability µ 0 , as shown in Sec. 5.As mentioned in the introduction, these two definitions differ from the definitions used in the metamaterial theory, but there is a direct link between them (see Sec. 5).
The average over the whole cell in Eq. ( 1) is related to the average over each region (1 and 2) by Here the subscript V (1,2) stands for the domain of averaging and, in general, ε eff is a tensor.In the same manner: On the other hand, from the macroscopic Maxwell equation we have the two standard relations, valid also for the mean fields: Substitution of Eqs. ( 2) -( 4) into Eq.( 1) results in the following relation: In order to obtain a form of ε eff that does not depend on the mean electric fields, it is necessary to explicitly determine, if possible, the relation between the mean fields: There are many different ways to determine this link, all based on the assumption that the volume of inclusion is so small that the fields are uniform, so that the link between the averaged values is the same as for the uniform fields.Two main approaches are usually used: 1) from the relation between the polarizability of the inclusion and the effective permittivity (Causius-Mossoti formula [8,9]) and 2) from the Green tensor singularities [10].Some less known methods use the boundary conditions in 1D case, or a matrix approach to solving Maxwell equations [11].In order that the effective media approach results in a linear effective permittivity, it is necessary that the tensor m Q does not depend on the electric field, so that In the case of highly symmetrical inclusions the Q m tensor is diagonal, so the above equation is simplified into Appendix A presents a Green tensor approach to obtain a tensor Q relating the field in the medium 1 to that in the medium 2 at the interface (see Eq. ( 61)).The resemblance of Eqs. ( 6) and (61) must not be misleading, because the electric fields in the two equations are generally quite different, as well the tensors that serves as a link.However, to combine the two approaches and derive expressions of the effective permittivity, we can suppose that the field in the two media is almost a constant, or in other words, well homogenized.This is equivalent to making the following strong assumption: Then, Eqs. ( 61) and (74) from Appendices A and B permit immediately to obtain the wellknown electric response of a spherical or cubic inclusion (or cavity) [1]: resulting in the well-known expression for the effective permittivity: For an infinitely long (in z direction) circular cylinder, by using Eqs.( 61) and (75), the link between the two fields is given as so that As far as we are able to find, in the case of cylindrical elliptical inclusions, an equation having Eq. ( 8) as a special case was first established by Nicorovici and McPhedran Eq. (39) of [12], and it corresponds to the conductivity relation of Galeener [13].
For a slab infinite in y and z directions, the link between the field values at the interface between V 1 and V 2 is determined from Eq. (76) and corresponds to the well-known boundary condition Equations ( 14) and ( 8) result in uniaxial effective permittivity.In direction parallel to the plane (TE polarization) it is given by the arithmetic mean value: In direction perpendicular to the plane (TM polarization) it is equal to the harmonic mean value: Equations ( 15) and ( 16) give the components of the well-known effective permittivity tensor of a 1D grating periodic in x direction in the limit of small period-to-wavelength ratio.
In the limit of much larger (in absolute value) permittivity (almost infinitely conducting metallic gratings) the last formula takes the form: As discussed after Eq. ( 6), the above formulas are obtained using a strong assumption m Q = Q .For simple geometries and for more complex ones (elliptical or rectangular cylinder) the relevant formulas are given in Appendix B. However, as shown in the next section, when the feature dimensions are comparable with the metal skin depth, this assumption is not valid and it is necessary to replace it by another one.

1D lamellar grating
In normal incidence and TM polarization the propagation constant of the fundamental mode that can propagate inside the grooves is proportional to eff ,xx ε . It is important to stress that Eq. ( 17) is valid in the limit of the assumption of Eq. ( 9), i.e. that the electric field in regions 1 and 2 is sufficiently homogenized.This occurs for two situations: (1) for dielectrics materials, in which the field is generally slowly varying, (2) for metals, when the skin depth L s is much greater than the metal width Φ, 2 2 Im(n )/ 1 πΦ λ << .In the latter case if |ε 2 | tends to infinity, the limit of eff ,xx ε is given by Eq. ( 17).
On the other hand, when the skin depth is much less than the metal width, in particular, when metals can be considered as perfectly conducting as in the microwave region, we can assume that E 2 = 0, as it is observed in wire polarizers, so that eff ,xx 1 The problem with these simple Eqs. ( 15)-( 17) is that they represent the static limit with 2 2 Im(n )/ 0 πΦ λ → , whereas in the visible or near IR the real structures are far away from this limit.As mentioned in the Introduction, one can find many attempts to overcome the limitation for 1D, 2D, 3D periodic structures, and for aperiodic distribution of inclusions.A detailed historical analysis can be found in [11].In case of 1D gratings, Rytov [14] proposed a correction for 1D grating, proportional to the second order of period-to-wavelength ratio (d/λ) 2 .Rytov obtained the following expression valid within second order of the small parameter: ( ) In TM polarization [14,15], present another closed-form expression: ( ) As we see further on, these two expressions work quite well when the wavelength is fixed and the period is reduced.However, they differ from the numerically obtained values even for dielectric gratings, when 2 2 Im(n )/ 0.5 πΦ λ > , which is quite natural, because the electric field varies significantly within one grating period.Moreover, for metallic gratings, when d is fixed and the wavelength increases while the metal dispersion is taken into account, the absolute value of ε 2 also increases and in many cases the expressions do not converge, as the wavelength tends to infinity.In the case of metallic gratings, first, the limitations of the validity of the expression (in TM polarization) with respect to d/λ ratio are much stronger.Moreover, they cannot take into account the case when the skin depth is smaller than the metal thickness.
Let us consider a lamellar metallic grating made of silver.Below the free-electron plasma frequency (above wavelength of 0.350 µm) the skin depths L s does not decrease with the wavelength and stays around 0.01 µm up to wavelength values exceeding 10 µm.16) in blue, Eq. ( 20) in red, and the correction by using Eq. ( 25) in black.
As a "good" example, we take the period equal to 30 nm and the metal lamella width w = 15 nm, which is comparable to the skin depth in the infrared.given by Eqs. ( 16), (20), and the numerical values obtained by the Fourier-modal method [16] and equal to the square of the eigenvalue of the mode inside the grating structure having the smallest imaginary part.A very good coincidence is observed between numerical results and the approximate formula, Eq. ( 20), even in the spectral domain where silver changes from a conducting metal (0.4 µm) to a lossy dielectric (0.3 µm).The mean harmonic values in Eq. ( 16) start to divert systematically from the values given by the other methods at longer wavelength, because w is comparable to L s .For the same reason, the results of Eq. ( 20) slightly differ from the numerical values and from the black line, that is made by introducing a correction to Eq. ( 16), described in the next paragraph, Eqs. ( 23) and (24).Equation ( 16) was obtained by applying the strong assumption of Eq. ( 9).If the period d and the width w are decreased 10 times more, its results coincide with the numerical values.However, it is clear that this assumption cannot work if the electric field is strongly varying within the grating period.The Fourier-modal method is rapidly converging and does not require more than +/− 10 Fourier harmonics.For covering partially this case, there is a limited number of possibilities.Concerning the Green's tensor approach, there are several options: 1.In addition to the singular part, it is possible to integrate the principal part, assuming small variation of the electric field E. Detailed analytical results for 2D to 3D cases can be found in [17].However, our numerical experiment showed that the corrections to the effective permittivity are quite large.In 1D case, Eq. ( 70) of Appendix B contains only a singular part for G xx , and the integration of G yy gives results that are of first order with respect to d/λ, a correction larger than necessary to fit the numerical results.The reason is that the hypothesis that E varies more slowly than the Green tensor is not valid.
2. To solve the integral equation with a full kernel (singular and principal value) with unknown E. However, it is more rigorous to use the Green tensor of the corresponding periodic structure, and the difficulties are the same.This leads to the rigorous integral method that already exists, and requires numerical solutions.
3. Intermediate hypothesis, that one can limit oneself to the singular part L, and to assume that the electromagnetic field varies in some known manner inside the medium 2. In what follows, we chose the third option by proposing a second, weaker assumption instead of Eq. ( 9).
The assumption is based on the physical fact that the field inside the metal decreases exponentially in direction perpendicular to the surface with a decay factor equal to 2 2 Im(n ) / π λ, and it vanishes beyond a skin depth L s distance, whereas it is much less rapidly varying inside the dielectric regions.In order to take this into account, we consider an approximation of the modal method in which the electromagnetic field in each medium is represented as a sum of modes, each one being an exact solution of the propagation equation.
Preserving only the fundamental mode, its x-dependence is given by the x-propagation so that at the boundaries 2 2 , 0 E (w / 2) E ≈   .If E 2,0 is much slowly varying than the cosine term, the mean field is given by: From Eqs. ( 6) and ( 22) we find that in this weaker assumption which, for lamellar grating in TM polarization becomes The black lines in Fig. 2 are obtained by using Eq. ( 8) together with Eq. ( 25).In the TE case In order to investigate the validity of the method for larger lamellae dimensions, Fig. 3 shows the same dependencies as in Fig. 2, but with grating dimensions 10 time larger, d = 0.3 µm and w = 0.15 µm, much more realistic than the previous ones.Here, the lamellae width is 15 times the skin depth in the infrared, so that we can expect that the field is not at all homogenized inside the metal, even for very large λ/d ratios.We can observe that the corrected approach accounts quite well not only for large wavelengths, but also when λ approaches d.Again, the corrected approach given by Eqs. ( 8), (25), and (26) works much better than the approximation (d/λ) 2 .The latter does not converge with increasing λ, because the absolute value of permittivity of Ag increases more rapidly than λ 1 , Fig. 6.As a result the right-hand side of Eqs. ( 19) and ( 20) tend to infinity as λ grows.
Several remarks are due: 1.When the lamellae width w is small, the correction factor q tends towards 1, i.e. to the strong assumption, and the static Eqs. ( 15) and ( 16) are valid.
2. For thicker conducting lamellae the correction factor q decreases, reflecting the fact that the mean field in region 2 becomes smaller than on the interface.For highly conducting gratings, q decreases as 3.One may ask why we have taken cosine dependence in Eq. ( 21).In case of normal incidence, electric field is symmetrical with respect to the origin, which is taken to be in the middle of region 2.An additional argument is that the antisymmetrical sine dependence gives a zero contribution to the mean field, when integrated from -w/2 to + w/2, so that we are using the symmetric cosine dependence.
Fig. 6.Spectral dependence of the real part of silver relative permittivity [18].
4. This renormalization approach cannot be applied to dielectric grating, because the tangent function diverges at odd multiple π/2, and vanishes at the even multiple of π/2.In the metallic case, they do not influence the results when the imaginary part of the refractive index is large.
5. Below λ = 0.35 µm, there appear other eigenvalues that have lower imaginary part than the one presented in Fig. 3, so that it is not reliable to consider that there is a single transfer channel inside the corrugated region.

2D gratings
In what follows, we present a study of gratings having 2D periodicity with a square unit cell and cylindrical inclusion, invariant in the vertical direction.We consider inclusions of square, rectangular, circular, and elliptical cross sections.In addition to Eq. ( 11) that works perfectly in the static limit, there are several other approaches that determine a lower and an upper limit of the effective permittivity.A detailed demonstration can be found in [19], where the authors consider the case of the mean electric conductivity.The lower bounds are obtained by considering the unit cells arranged at first in series in one direction, summing the resistivity, and then summing the conductivity (inverted resistivity) in the other direction.The upper bounds are obtained by considering the arrangement at first in parallel by summing the conductivity in one direction, then summing the reciprocal (i.e., the resistivity) in the other direction, and then inverting the result to obtain the conductivity.These two bounds can be directly applied for the permittivity, although the inverse of the permittivity has no direct physical sense.
It is interesting to note that the two bounds obtained in [19] can be derived from the Fourier modal method when applied to 2D gratings with rectangular inclusions.For 2D gratings, there can be many equivalent Fourier representations of the permittivity function, all obeying the Fourier factorization rules.Among them the formulation in [16] gives the lower bound and that in [20] gives the upper bound, both in the form of the 00-th element of the permittivity function.Therefore, Here f x and f y are the filling ratios in x-and y-directions, respectively.Figure 7 presents the dependence of effective permittivity on the wavelength for an array of silver circular cylinders, when their dimensions are much smaller than the wavelength and the skin depths.The side of the unit cell is d = 3 nm and the cylinder diameter is 2R = 1.5 nm., with a filling factor f = πR 2 /d 2 .As can be observed in Fig. 7(a), the static limits is included within the two bounds, and coincide with the numerical values, even in the resonant domain around the free-electron plasma frequency, where the upper and the lower bounds are inverted, Fig. 7(b).Contrary to the 1D case, the Fourier modal method in the 2D case has slower convergence rate and it is necessary to use not less than 25 Fourier harmonics in each direction, i.e. totally 51x51 harmonics, which requires much longer computation time (several minutes per point).
As can be expected from the 1D case, if the feature dimensions are larger than the skin depth, the static approach is not valid, unless perfect conductivity is assumed.And indeed, Fig. 8 shows that the numerically obtained values are lower than both the bounds and the static predictions, when real metal is considered.On the contrary, if a correction to Eq. ( 8) is introduced in the same manner as in Eqs. ( 23) and (24), it is possible to quite well match the numerical results, as observed in Fig. 8.
In cylindrical geometry, the natural assumption of the field behavior inside the metallic cylinder is that it is proportional to the Bessel function of the 1st kind and 0th order J 0 .With an almost purely imaginary argument for highly conducting metals, J 0 is equal to the modified Bessel function I 0 , and decreases rapidly in the metal depth.The coefficient of proportionality has to be taken, as in Eq. ( 21), to be equal to the value of J 0 at the cylinder surface, ρ = R: Fig. 8. Same as in Fig. 7, but with d = 300 nm and R = 75 nm.The black curve is made by using the correction in Eq. (29).
As in 1D case, assuming that E 2,0 varies in ρ much weaker than J 0 , the mean electric field over the cylinder surface is then given by: As in the 1D case, the correction factor tends to 1 when R becomes small compared to λ/|n 2 |, i.e., smaller than the skin depth.The black curve in Fig. 8 is from Eq. ( 8) using this correction q to Eq. ( 24).Let us remind that the singularity L is diagonal and isotropic in the x-y plane, xx = L yy = 1/2, L zz = 0, Eq. ( 75).In the region below 0.4 µm, there appear a multitude of eigenvalues with small imaginary part, see Fig. 9, that allow no definitive answer which approach is better.In addition, the existence of these eigenvalues show in the same manner as in the 1D case, that when λ approaches the period, any approach that lays on a single channel of transmission by the structure cannot succeed.
The singularity tensor for cylinders with square section is the same as for circular cylinders, so that Eq. ( 8) can be applied directly, by taking into account the difference in the filling factor f = w 2 /d 2 .Figure 10 presents numerical and different analytical results for square and circular cylinders as a function of the filling factor f at wavelength of 2 µm.As can be observed, for very small period, the numerical results match perfectly the static formula, and lies between the two bounds.When the dimensions become larger than the skin depth, numerical results differ significantly from the static ones.On the contrary, the modified formula by taking into account Eq. ( 29) gives results closer to the numerical values, the difference growing with the filling factors.) on the filling factor f for an array of circular or square cylinders.Dots present the numerical results (red squares with d = 3 nm, black for squares, and green for circles with = 300 nm).Red curve Eq. ( 13), lower and upper limits, Eq. ( 27) in blue and rose, and the correction from Eq. (29) in black.
The general tendency observed in the examples shows that in the non-homogenized cases, the effective permittivity has lower real part than in the static case, because in the former, the electric field penetrates only inside the skin depth in the metal and is more localized in the surrounding dielectric (air) with n = 1.

Cylinders with elliptic and rectangular cross-sections
The singularities for elliptic profiles can also be obtained analytically (Appendix B), and they introduce anisotropy in the x-y plane.In order to obtain a correction similar to Eq. (29) in the static equation for the effective permittivity, we have several requirements: (1) the correction function must be equal to 1 over the elliptical surface of the inclusion and ( 2) it has to vary in the direction normal to this surface.There exist elliptic functions with these properties that are solutions of propagation equation in elliptical coordinates, Mathieu functions.However, we were not able to find their analytical integration, and we shall use an approximate approach.A natural generalization of Eq. ( 29) that provides analytical formulas is given by the Bessel function, with an argument that is constant over the ellipse: where R is a free parameter.We chose it so that the surface of the equivalent circle is equal to the elliptic surface: Equation (30) satisfies the two requirements, it is constant on the ellipse and its change is presented by the gradient of the function, which is perpendicular to the ellipse.The integration of Eq. (30) over the ellipse surface is straightforward: which is exactly the same as for a circle, only that the value of R is given by (32).We can use the same approach for a rectangular cross-section, with the singularities given in Appendix B, Eq. ( 81).The comparison between the numerical data and approximate formulas gives results, similar to those for circular cylinders, as observed in Fig. 10.

Effective surface impedance and the reflection of the homogenized layer
As remarked by Rytov at the end of his paper, while in the depth of the grating region it is usually sufficient to consider a single transmission channel (having the smallest imaginary part of the constant of propagation in z-direction) between the cladding and the substrate, on their boundaries even the fast decaying field components can play an important role.And indeed, Fig. 11 shows that in the case of a 1D lamellar grating, the effective permittivity of the fundamental transmission converges much more rapidly when increasing the number of Fourier harmonics in the Fourier modal method than the reflectivity of the grating.
The precision of 1% for ε eff is reached by taking into account no more than ± 3 Fourier harmonics, whereas for the reflectivity it is necessary to use at least ± 10 harmonics.
If the grating region is considered as an equivalent plane anisotropic but homogeneous layer, the reflection coefficients at each of its plane boundaries are given by the Fresnel In normal incidence, the reflection between the cladding and the homogeneous layer is simply given by, for x-polarized electric field where n s is the ratio between the magnetic and electric field at the boundary between the cladding and the grating region (i.e. the normalized inverse impedance), and n cl is the cladding refractive index.The reflection coefficient of the layer with thickness h is then given by Airy formula: The most direct choice is to take the reciprocal relative impedance given simply by the normalized propagation constant in z-direction: In the case of completely homogenized field, this formula gives the same results as the numerical results for the periodic structure, as shown in red in Figs.12(a) and 13(a).However, in case of larger period and lamellae width, the numerical values of the reflectivity are almost three times larger than when using the "homogenized" values of n s .Moreover, as observed in the results obtained all over the paper in the non-homogenized cases (feature dimensions larger than the skin depth), the real part of the effective permittivity is smaller than in the homogenized case, thus the reflection of the equivalent layer will be even smaller than the numerical values.The reason is that the approach in Sec.2 and 3 requires the continuity of the tangential E and normal D components of the field, which is not the case close to the cladding/grating (and substrate/grating) interfaces.On the sharp edges E and D are not only discontinuous, but also diverging.In order to obtain satisfactory approximation, we use the technique developed in the exact modal method [21,22].
Let us consider the interface the semi-infinite homogeneous cladding region and a semi-infinite grating region with the interface at z = 0.In TM polarization, the magnetic field in the cladding is represented as a sum of plane waves, and in the grating structure as a sum of modes, each one exact solution of Maxwell equations.16) and (36).d = 300 nm and w = 150 nm, in black, numerical results, in blue, modal method using Eq.(42) for n s .(b) dependence of the first maximum of the reflectivity on the lamellae width w for d = 300 nm (n 2 = 0.129 + i 6.83 [18]).
For our purpose, we truncate this representation to single incident and reflected plane waves in the cladding and a single downward propagating mode in the grating: where In the modal method different projection methods can be used to solve the boundary matching equations.Here we use one of the hybrid methods [21][22][23].The equation resulting from the continuity of H at z = 0 is projected on the basis of the plane waves in the cladding: whereas the continuity of the field normal (in z) derivative is projected on the basis of the modal functions (x) ψ : ( ) Similarly to Sec. 2, in region 1 filled with low-index dielectric, the field can be considered constant, so that where the single and double primes stand for the real and imaginary parts, respectively.These integrals represent approximations of the more complex overlap integrals presented in [23].
The blue line in Fig. 12 presents the results of the reflectivity obtained from Eqs. ( 34) and (35) with n s given from Eq. (42).As can be observed, the reflectivity values match the numerical results much better than by using Eq.(36) for the entire region of lamellae width.In addition, when w becomes much larger than the skin depths, n s from Eq. ( 42) tends towards the static value of Eq. (36).
In case of 2D periodicity the tendency is the same as in the 1D case, as can be seen in Fig. 13: 1.For very small feature size, the static approach is correct for both ε eff and n s from Eq. (36).
2. For larger features the reflectivity is higher than the static limit, and the period in h is longer, because eff,xx Re(n ) is smaller when using the correction in Eq. (29).
In order to obtain a relevant correction, we apply the same approach as for the 1D case.Let us consider an incident plane wave with its magnetic vector oscillating in x-direction.The magnetic field in the cladding is represented as a sum of plane waves: In the grating region we assume that the field can be factorized Using the same procedure as for the 1D case, we obtain the following relations: that enables us to derive the same relation as Eq. ( 42).As in Eq. ( 33) the integral of ψ can be evaluated analytically: The second integral can be explicitly presented in the following form, provided that 2 2 Re n Im n << :

Metamaterial parameters and effective medium permittivity and permeability
As explained in the introduction, in the paper we follow the classical definitions of effective permittivity, Eq. ( 1), and effective permeability, given by a similar relation, definitions essentially used during the 20th century, not only for the permittivity, but also for the conductivity [6,7].For non-magnetic materials, thus defined effective permeability is equal to the vacuum permeability µ 0 .A simple demonstration follows the logic of Eqs. ( 1)-( 5).The average permeability gives the link between the averaged magnetic field and induction: For non-magnetic materials the local permeability is everywhere equal to the vacuum permeability, and because of the linearity of the integration: whatever the special distributions of the fields, so that it directly follows that It is necessary to stress out that the definition of effective permittivity and permeability used in this paper represent the link between the averaged fields and differ from the definitions used in most of the papers devoted to metamaterials [24,25], although there is a direct link between them.The metamaterial effective permittivity M ε and permeability M μ are obtained from the fundamental Bloch mode propagation constant n M and the reflection coefficient r cl on the interface between the metamaterial and the homogeneous layer, defined in Sec. 4. In the case of n cl = 1, the link is given by: In most of the metamaterial studies the values of n M and r cl are calculated using some numerical method, except for ref [23].This allows to determine M ε and M μ using Eq. ( 53).
On the contrary, in Sec. 2 and 3 we were able to introduce an analytical correction to the average permittivity obtained from the homogenization approach, summarized in Sec. 1.This correction permits to directly obtain the fundamental Bloch mode propagation constant: for non-magnetic materials.Subsequently, the correction to r cl coming from Eqs. ( 34) and (42) enables us to analytically obtain its values, so that the metamaterial permittivity and permeability can be determined from where n s is the normalized inverse impedance.In the case of (almost) completely homogenized fields and non-magnetic materials, In the cases of non-homogenized the fundamental mode propagation constant is different from the optical index participating in the surface impedance, eff s n n ≠ , and even non-magnetic materials can have magnetic response with M 1 μ ≠ , although the average permeability as defined in Eq. ( 50) is equal to one.Figure 14 presents the spectral dependence of the real parts of metamaterial permittivity and permeability for the case of 1D lamellar grating having two different periods, 3 nm with homogenized fields, and 300 nm.As observed, in the small-feature case the permeability calculated from Eq. ( 55) is equal to one, whereas for large feature-size the system shows magnetic properties.
Another example concerns an Ag grating suspended in air with a period d = 0.15 µm and channel width of 0.0256 µm (w = 0.1244 µm) that shows several anomalies in reflection in the visible, as observed in Fig. 15(a).The numerical results (in black) are compared with the reflectivity obtained by using the static values of n eff and n s (red curve), and with their corrected values from Sec. 2 and 4 in blue.The corrected results match quite well the numerical values, even taking into account that λ/d ratio is relatively small (from 2 to 4).
If the period and the lamella width are reduced 100 times, the static limit reflectivity match the numerical results (red circles), because the dimensions are much larger than the skin depth and the fields are completely homogenized.55).They vary strongly within the spectral interval and, in particular, the structure shows well pronounced magnetic properties.It is necessary to stress that the values of n eff and n s are not found numerically, but using the analytical corrections to the static limit from Sec.2 and 4.
On the contrary, the small-feature structure (d = 1.5 nm) is non-magnetic, as shown with red circles in Fig. 15(b).

Conclusions
The phenomenological physical hypothesis for the electromagnetic field behavior inside finite conducting metal features of 1D and 2D gratings having invariant vertical geometry serves to obtain both effective permittivity and surface impedance in the cases when electromagnetic field is not sufficiently homogenized to behave according to the static approximation.In particular, is the case in the visible and near IR if the metallic feature size is comparable with the metal skin depth.
The analytical expressions for the correction factor are also valid in the two limiting wellknown cases: 1.In the static limit, when the grating dimensions are so small with respect to the skin depth that the electromagnetic field can be considered constant (correction factor q tends to 1).
2. For perfectly conducting metals (q tends towards 0).They also permit us to determine the metamaterial permittivity and permeability, as defined in Sec. 5.In particular, our approach confirms that non-magnetic gratings can show magnetic metamaterial behavior in the case of non-homogenized fields, even if the relative effective permeability, defined as the ratio between the average induction and the average magnetic field is equal to one.

Appendix A. Green's tensor relation between the electric fields in V 1 and V 2
Let us consider at first the Green's function approach to light scattering by an object having permittivity ε 2 embedded in a medium with permittivity ε 1 .
The total field p E (r)   is the sum of the field solution of the unperturbed problem unp E (r)   (obtained if V 2 =0) and the field scattered by the inclusion (volume V 2 ): where the bold one stands for the three-dimensional unit tensor, and the electric dyadic Green's tensor is the solution of The scattered field depends on the total field all over the scattering object through the electric dyadic Green's tensor G.This tensor has a singular part L and a principal-value part P vG : where (r r ) ′ δ −    is the Dirac function.In the first-order approximation, the strongest contribution is the auto-scattering part expressed by the singular part of the Green's function: Wherefrom at a point S r  on the boundary S between the two media: The trace of L must be equal to one: Trace( ) 1.

= L
The singular part is quite often called depolarization factor and it has a simple form and depends on the form of the object.A detailed analysis and many examples can be found in [10].In Appendix B, we present a short review of Green tensor singularities.

Appendix B. Green's tensor and its singularities
For the sake of completeness, here we summarize some of the properties of Green's tensor singularities that can be found dispersed among many different publications.For the electric field, the Green's tensor is a solution of the singular propagation equation: Equation ( 64) can be further developed: On the other hand, the scalar Green's function satisfies the equation: where the derivatives are taken with respect to the observation coordinate r  .From this equation it becomes clear that even in the case of no singularity of the Green's function, the tensor can have singularities.For example, in the 1D case: The singularities of G are treated by integrating them in a volume of the inclusion V 2 , excluding a small volume V 0 around the singular point, with V 0 tending to zero.Applying the theorem of Ostrogradski-Gauss, the volume integral is transformed to a surface integral.In the electrostatic limit the 3D Green's function is where n is the unit vector normal to the surface.
In the 2D case and the source term positioned in the origin, the static potential is given by lnρ/2π, and the integration for L is carried over the contour of inclusion C 2 : For a sphere and a cube, the symmetry of the system results in a symmetrical depolarization factor: 1/ 3 0 0 = 0 1 / 3 0 .
For a long cylindrical inclusion with circular or square cross-section (see Secs.B.1 and B.2) having axis along z-direction: 1/2 0 0 = 0 1/2 0 .0 0 0 For a 1D grating, periodical in x-direction: 1 0 0 = 0 0 0 .0 0 0 Finally, for an infinitely thin circular disk the singularity of the Green's function tensor takes the form: And the link between the total and the incident field takes the following form: i.e. the electric field tangential components and the displacement normal component are continuous on the disk surface, as expected.
A demonstration for rectangular and elliptical cross-sections is given below:

Fig. 1 .
Fig. 1.Schematic presentation of the elementary cell of a 1D (a) or 2D (b) and (c) periodical system with notations.

Figures 4
Figures 4 and 5 present the dependence of the effective permittivity in TE and TM polarization for λ/d = 3 and the different approaches as a function of the filing factor f = w/d.Again, the corrected approach given by Eqs.(8),(25), and (26) works much better than the approximation (d/λ) 2 .The latter does not converge with increasing λ, because the absolute value of permittivity of Ag increases more rapidly than λ 1 , Fig.6.As a result the right-hand side of Eqs.(19) and (20) tend to infinity as λ grows.Several remarks are due:

Fig. 11 .
Fig. 11.Convergence of the relative error of the real part of ε eff and of the reflectivity as a function of the number of Fourier harmonics [-N, N] for a lamellar diffraction grating made of silver with d = 0.3 µm, w = 0.15 µm suspended in air in normal incidence and TM polarization, λ = 1 µm, h = 0.25 µm.cl sub 0 eff cl sub 0 eff r r exp(2ik n h) r , 1 r r exp(2ik n h) + = + (35)where r sub is the reflection coefficient at the substrate boundary and eff

Fig. 12 .
Fig. 12. Reflection by a 1D lamellar silver grating suspended in air in normal incidence and TM polarization at λ = 1 µm.(a) dependence on the grating thickness h.In red, d = 3 nm, w = 1.5 nm, with points, numerical results, and with a line, Eq. (35) using Eqs.(16) and (36).d = 300 nm and w = 150 nm, in black, numerical results, in blue, modal method using Eq.(42) for n s .(b) dependence of the first maximum of the reflectivity on the lamellae width w for d = 300 nm (n 2 = 0.129 + i 6.83[18]).
and (40) enable us to obtain the formula for n s :

Fig. 13 .
Fig. 13.(a) Reflectivity of a 2D array of silver circular cylinders (Fig. 1(c)) in air as a function of the grating height h.Normal incidence, λ = 1 µm.Red points, numerical values with d = 3 nm, R = 1 nm, red curve, the static limit.Black curve, numerical results for d = 300 nm, R = 100 nm, blue line, correction from Eqs. (42), (48), and (49).(b) the dependence of the maximum of the reflectivity on the cylinder diameter for d = 300 nm.
the previous sections we used the relation eff eff ,xx n = ε .

Fig. 14 .
Fig. 14.Spectral dependence of the real parts of the metamaterial permittivity ε M and permeability µ M , calculated using the values of n eff and n s from Sec. 2 and 4 for a lamellar Ag grating with 50% filling ratio, and having two different periods, (a) d = 3 nm in red and (b) d = 300 nm in black.
discussed in connection with Figs. 2, 7, and the red lines in Figs. 12 and 13 and thus M M e f f , x x , 1. ε = ε μ =

Fig. 15 .
Fig. 15.1D lamellar Ag grating in air with d = 0.15 µm, h = 0.24 µm, and w = 0.1244 µm in TM polarization.(a) Reflectivity as a function of the wavelength, numerical results in black, modal correction in blue, the static results, red line.The numerical results for d = 1.5 nm and w = 1.244 nm is presented with red circles.(b) The spectral dependence of the metamaterial parameters ε M and µ M , calculated using Eq.(55) with the modal correction for n eff and n s from Sec. 2 and 4.

Figure 15 (
Figure15(b)  presents the dependence of the metamaterial permittivity and permeability obtained using Eq.(55).They vary strongly within the spectral interval and, in particular, the structure shows well pronounced magnetic properties.It is necessary to stress that the values of n eff and n s are not found numerically, but using the analytical corrections to the static limit from Sec.2 and 4.On the contrary, the small-feature structure (d = 1.5 nm) is non-magnetic, as shown with red circles in Fig.15(b).
of this equation and taking into account that div(rot) = 0

0 1 /
4π − r r   , i.e., it is necessary to integrate: of the inclusion V 2 .Assuming that the source term is in the origin: