1. 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.

2. 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).

 

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

Download Full Size | PDF

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 ${\epsilon }_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, ${\epsilon }_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₂/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:

The average over the whole cell in Eq. (1) is related to the average over each region (1 and 2) by

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:

For a slab infinite in y and z directions, the link between the field values at the interface between V₁ and V₂ is determined from Eq. (76) and corresponds to the well-known boundary condition

In the limit of much larger (in absolute value) permittivity (almost infinitely conducting metallic gratings) the last formula takes the form:

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.

3. Numerical examples and counter-examples, modified approach

3.1. 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 $\sqrt{{\epsilon }_{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\pi \Phi \mathrm{Im}(n_2)/\lambda <<1$. In the latter case if |ε₂| tends to infinity, the limit of ${\epsilon }_{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₂ = 0, as it is observed in wire polarizers, so that

The problem with these simple Eqs. (15)-(17) is that they represent the static limit with $2\pi \Phi \mathrm{Im}(n_2)/\lambda \to 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/λ)². Rytov obtained the following expression valid within second order of the small parameter:

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.

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. Figure 2 presents the real part of ${\epsilon }_{eff,xx}$ 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.

 

Fig. 2 Spectral dependence of Re(${\epsilon }_{eff,xx}$) for a 1D lamellar grating (the medium 2 in Fig. 1(a) is silver. Medium 1 is air, d = 30 nm, w = 15 nm). Long (a) and short (b) wavelength regions. Numerical values (Fourier modal method) are presented with dots, harmonic mean from Eq. (16) in blue, Eq. (20) in red, and the correction by using Eq. (25) in black.

Download Full Size | PDF

For covering partially this case, there is a limited number of possibilities. Concerning the Green’s tensor approach, there are several options:

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\pi \mathrm{Im}(n_2)/\lambda$, 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 constant $k_{2,x}=\sqrt{k_2^2-k_0^2{\epsilon }_{eff}}$ with k₀ is the vacuum wavenumber. For highly conducting gratings $k_{2,x}\approx k_2$, and fixing the coordinate origin in the middle of medium 2, we assume that

 

Fig. 3 Same as in Fig. 2, but d = 300 nm and w = 150 nm.

Download Full Size | PDF

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/λ)². The latter does not converge with increasing λ, because the absolute value of permittivity of Ag increases more rapidly than λ¹, Fig. 6. As a result the right-hand side of Eqs. (19) and (20) tend to infinity as λ grows.

 

Fig. 4 Re(${\epsilon }_{eff,yy}$) for the case presented in Fig. 3 as a function of the filling factor f = w/d. λ = 1µm.

Download Full Size | PDF

 

Fig. 5 Same as in Fig. 4 but for TM polarization.

Download Full Size | PDF

 

Fig. 6 Spectral dependence of the real part of silver relative permittivity [18].

Download Full Size | PDF

Several remarks are due:

3.2. 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,

Figure 7 presents the dependence of the 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²/d². 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).

 

Fig. 7 Spectral dependence of Re(${\epsilon }_{eff,xx}$) = Re(${\epsilon }_{eff,yy}$) for an array of silver circular cylinders in air, Fig. 1(c). Square cell with d = 3 nm and R = 0.75 nm. Numeric values presented with dots, static limit, Eq. (13) in red, lower and upper limits, Eq. (27) in blue and rose. Long (a) and short (b) wavelength domains.

Download Full Size | PDF

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.

 

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).

Download Full Size | PDF

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₀. With an almost purely imaginary argument for highly conducting metals, J₀ is equal to the modified Bessel function I₀, 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₀ at the cylinder surface, ρ = R:

As in 1D case, assuming that E_2,0 varies in ρ much weaker than J₀, the mean electric field over the cylinder surface is then given by:

 

Fig. 9 Same as in Fig. 8 but for shorter wavelengths.

Download Full Size | PDF

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²/d². 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.

 

Fig. 10 The dependence of Re(${\epsilon }_{eff,xx,yy}$) 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 d = 300 nm). Red curve Eq. (13), lower and upper limits, Eq. (27) in blue and rose, and the correction from Eq. (29) in black.

Download Full Size | PDF

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.

3.3. 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:

The integration of Eq. (30) over the ellipse surface is straightforward:

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.

4. 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.

 

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.

Download Full Size | PDF

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 formulas. In normal incidence, the reflection between the cladding and the homogeneous layer is simply given by, for x-polarized electric field

The most direct choice is to take the reciprocal relative impedance given simply by the normalized propagation constant in z-direction:

 

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₂ = 0.129 + i 6.83 [18]).

Download Full Size | PDF

 

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.

Download Full Size | PDF

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 between 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.

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:

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:

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:

5. 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 µ₀. A simple demonstration follows the logic of Eqs. (1)-(5). The average permeability gives the link between the averaged magnetic field and induction:

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 ${\epsilon }_M$ and permeability ${\mu }_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:

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:

In the cases of non-homogenized fields the fundamental mode propagation constant is different from the optical index participating in the surface impedance, $n_{eff}\ne n_s$, and even non-magnetic materials can have magnetic response with ${\mu }_{{}_M}\ne 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.

 

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.

Download Full Size | PDF

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).

 

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.

Download Full Size | PDF

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.

Figure 15(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).

6. 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, this 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 well-known cases:

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₁ and V₂

Let us consider at first the Green’s function approach to light scattering by an object having permittivity ε₂ embedded in a medium with permittivity ε₁.

The total field ${\overrightarrow{E}}_p(\overrightarrow{r})$ is the sum of the field solution of the unperturbed problem ${\overrightarrow{E}}_{unp}(\overrightarrow{r})$ (obtained if V₂=0) and the field scattered by the inclusion (volume V₂):

(57)$${\overrightarrow{E}}_p(\overrightarrow{r})={\overrightarrow{E}}_{unp}(\overrightarrow{r})+k_0^2{\epsilon }_1{\displaystyle \underset{V_2}{\int }G(\overrightarrow{r}-\overrightarrow{r^{\prime }})\left(\frac{{\epsilon }_2}{{\epsilon }_1}-1\right){\overrightarrow{E}}_p(\overrightarrow{r^{\prime }})d\overrightarrow{r^{\prime }}},$$

where the bold one stands for the three-dimensional unit tensor, and the electric dyadic Green’s tensor is the solution of $rotrotG-k_0^2{\epsilon }_{1}G=1\delta (r-r^{\prime }).$

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:

(58)$$G(\overrightarrow{r}-\overrightarrow{r^{\prime }})=P_{vG}(\overrightarrow{r}-\overrightarrow{r^{\prime }})-\frac{1}{k_0^2{\epsilon }_1}L\delta (\overrightarrow{r}-\overrightarrow{r^{\prime }}),$$

where $\delta (\overrightarrow{r}-\overrightarrow{r^{\prime }})$ 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:

(59)$${\overrightarrow{E}}_p(\overrightarrow{r})\equiv {\overrightarrow{E}}_2(\overrightarrow{r})\approx {\overrightarrow{E}}_{unp}((\overrightarrow{r})-L\left(\frac{{\epsilon }_2}{{\epsilon }_1}-1\right){\overrightarrow{E}}_p((\overrightarrow{r}),in{V}_2$$

(60)$${\overrightarrow{E}}_p(\overrightarrow{r})\equiv {\overrightarrow{E}}_1(\overrightarrow{r})\approx {\overrightarrow{E}}_{unp}((\overrightarrow{r}),in{V}_1.$$

Wherefrom at a point $\overrightarrow{r_S}$ on the boundary S between the two media:

(61)$${\overrightarrow{E}}_2(\overrightarrow{r_S})\approx Q{\overrightarrow{E}}_1(\overrightarrow{r_S})$$

with

(62)$$Q={\left[1+L\left(\frac{{\epsilon }_2}{{\epsilon }_1}-1\right)\right]}^{-1}.$$

The trace of L must be equal to one:

(63)$$Trace(L)=1.$$

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:

(64)$$rotrotG-k^{2}G=1\delta (\overrightarrow{r}-{\overrightarrow{r}}^{\prime }).$$

By taking divergence of this equation and taking into account that div(rot) = 0, we obtain

(65)$$-k^2\nabla \cdot G=1\nabla \delta \Rightarrow \nabla \nabla \cdot G=-\frac{1}{k^2}\nabla \nabla \delta .$$

Equation (64) can be further developed:

(66)$$\nabla \nabla \cdot G-\Delta G-k^{2}G=1\delta (r-r^{\prime })$$

so that

(67)$$-\frac{1}{k^2}\nabla \nabla \delta -\Delta G-k^{2}G=1\delta \iff \Delta G+k^{2}G=-\left(1+\frac{1}{k^2}\nabla \nabla \right)\delta .$$

On the other hand, the scalar Green’s function satisfies the equation:

(68)$$\Delta g+k^{2}g=-\delta (\overrightarrow{r}-{\overrightarrow{r}}^{\prime })$$

and the relation between the Green’s function and tensor becomes

(69)$$G(\overrightarrow{r}-{\overrightarrow{r}}_0)=\left(1+\frac{1}{k^2}\nabla \nabla \right)g(\overrightarrow{r}-{\overrightarrow{r}}_0),$$

where the derivatives are taken with respect to the observation coordinate $\overrightarrow{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:

(70)$$\begin{array}{l}g(x-x_0)=\frac{i}{2k}\exp \left(ik\left|x-x_0\right|\right)\\ \frac{d^2}{d{z}^2}g(x-x_0)=-\delta (x-x_0)-k^{2}g(x-x_0)\\ G(x-x_0)=\left(\begin{array}{ccc}-\frac{1}{k_0^2{\epsilon }_1}\delta (x-x_0)& 0& 0\\ 0& g(x-x_0)& 0\\ 0& 0& g(x-x_0)\end{array}\right)\end{array}$$

The singularities of G are treated by integrating them in a volume of the inclusion V₂, excluding a small volume V₀ around the singular point, with V₀ 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 $1/4\pi \left|\overrightarrow{r}-{\overrightarrow{r}}_0\right|$, i.e., it is necessary to integrate:

(71)$$\nabla \nabla \frac{1}{4\pi \left|\overrightarrow{r}-{\overrightarrow{r}}_0\right|}$$

over the volume of the inclusion V₂. Assuming that the source term is in the origin:

(72)$$L=-{\displaystyle \underset{V_2}{\int }\nabla \nabla \frac{1}{4\pi \left|\overrightarrow{r}-{\overrightarrow{r}}_0\right|}}\to -{\displaystyle \underset{S_2}{\int }\widehat{n}\nabla \frac{1}{4\pi \left|\overrightarrow{r}-{\overrightarrow{r}}_0\right|}}dS={\displaystyle \underset{S_2}{\int }\widehat{n}\frac{\overrightarrow{r}-{\overrightarrow{r}}_0}{4\pi {\left|\overrightarrow{r}-{\overrightarrow{r}}_0\right|}^3}}dS,$$

where $\widehat{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₂:

(73)$$L={\displaystyle \underset{C_2}{\int }\stackrel{⌢}{n}\frac{\overrightarrow{\rho }}{2\pi {\rho }^2}dc}.$$

For a sphere and a cube, the symmetry of the system results in a symmetrical depolarization factor:

(74)$$L=\left(\begin{array}{ccc}1/3& 0& 0\\ 0& 1/3& 0\\ 0& 0& 1/3\end{array}\right).$$

For a long cylindrical inclusion with circular or square cross-section (see Secs.B.1 and B.2) having axis along z-direction:

(75)$$L=\left(\begin{array}{ccc}1/2& 0& 0\\ 0& 1/2& 0\\ 0& 0& 0\end{array}\right).$$

For a 1D grating, periodical in x-direction:

(76)$$L=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right).$$

Finally, for an infinitely thin circular disk the singularity of the Green’s function tensor takes the form:

(77)$$L=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right).$$

And the link between the total and the incident field takes the following form:

(78)$${\overrightarrow{E}}_2(\overrightarrow{r})\approx \left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& {\epsilon }_1/{\epsilon }_2\end{array}\right){\overrightarrow{E}}_1(\overrightarrow{r}).$$

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:

Infinite cylinder in z direction

B.1. Rectangular (square) inclusion

According to Fig. 16,

 

Fig. 16 Schematical view of a rectangular cross-section with notations.

Download Full Size | PDF

(79)$$on{c}_1:\{\begin{array}{l}x=\frac{a}{2}\\ \widehat{n}=\widehat{x}\\ \rho =x\widehat{x}+y\widehat{y}\\ d{c}_1=dy\end{array}$$

(80)$${\text{L}}_{\text{xx}}=\frac{2}{2\pi }{\displaystyle \underset{{\text{c}}_1}{\int }\frac{\text{xdy}}{{\text{x}}^2+{\text{y}}^2}}.$$

Using the substitution $y=\frac{a}{2}tg\phi$

(81)$$L_{xx}=\frac{1}{2\pi }{\displaystyle \underset{c_1}{\int }\frac{\frac{a^2}{4{\cos }^2\phi }d\phi }{\frac{a^2}{4}+\frac{a^2}{4}t{g}^2\phi }}=\frac{1}{\pi }{\displaystyle \underset{c_1}{\int }d\phi }=\frac{2}{\pi }arctg\frac{b}{a}$$

Similarly $L_{yy}=\frac{2}{\pi }arctg\frac{a}{b}$.

For squares b = a and $L_{xx}=L_{yy}=\frac{1}{2}$

B.2. Elliptical cross-section

If the semi -axes are equal to a and b, in parametric representation

(82)$$\begin{array}{l}\overrightarrow{\rho }=a\cos \phi \widehat{x}+b\sin \phi \widehat{y}\\ \widehat{n}=grad\rho /\left|grad\rho \right|\\ dc=\left|\frac{d\rho }{d\phi }\right|\end{array}$$

(83)$$\begin{array}{l}{L}_{xx}=\frac{1}{2\pi }{\displaystyle \underset{c}{\int }\frac{ab{\cos }^2\phi d\phi }{a^2{\cos }^2\phi +b^2{\sin }^2\phi }}=\frac{b}{a+b}\\ L_{yy}=\frac{1}{2\pi }{\displaystyle \underset{c}{\int }\frac{ab{\sin }^2\phi d\phi }{a^2{\cos }^2\phi +b^2{\sin }^2\phi }}=\frac{a}{a+b}\end{array}$$

Same as for squares, a circle is characterized by $L_{xx}=L_{yy}=\frac{1}{2}$.

Acknowledgment

Lifeng Li acknowledges the financial support of Aix-Marseille University for his recent one-month work visit to Institute Fresnel.

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]  

11. P. Lalanne and D. Lemercier-Lalanne, “On the effective medium theory of subwavelength periodic structures,” J. Mod. Opt. 43(10), 2063–2085 (1996). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]  

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). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]