Interferometric description of optical metamaterials

We introduce a simple theoretical model that describes the interaction of light with optical metamaterials in terms of interfering optical plane waves. In this model, a metamaterial is considered to consist of planar arrays of densely packed nanoparticles. In the analysis, each such array reduces to an infinitely thin homogeneous sheet. The transmission and reflection coefficients of this sheet are found to be equal to those of an isolated nanoparticle array and, therefore, they are easy to evaluate numerically for arbitrary shapes and arrangements of the particles. The presented theory enables fast calculation of electromagnetic fields interacting with a metamaterial slab of an arbitrary size, which, for example, can be used to retrieve the effective refractive index and wave impedance in the material. It is also shown to accurately describe optically anisotropic metamaterials that in addition exhibit strong spatial dispersion, such as bifacial metamaterials.

Optical metamaterials are man-made materials composed of densely packed subwavelengthsize nanoparticles appearing like artificial atoms to light. While the optical response of each individual nanoparticle can be revealed using, e.g. the electromagnetic multipole expansion [1], the description of the macroscopic optical response of a real three-dimensional metamaterial still remains a challenge. This description is complicated by non-trivial interactions between the nanoparticles, including evanescent-wave coupling between them. In this work however, we propose a way to calculate the transmission and reflection characteristics of metamaterial slabs, including anisotropic and spatially dispersive ones, without resorting to evanescent waves.
So far, both numerical [2][3][4][5] and experimental [6][7][8][9] techniques have been used to obtain the transmission and reflection coefficients for metamaterial slabs with a rather limited number of nanostructured layers. The coefficients obtained for such slabs do not necessarily describe the properties of a bulk metamaterial. One approach to characterize the material would be to successively increase the number of layers and see if the optical characteristics converge [3,4,9]. However, the understanding of the real physics that determines the final transmission properties is lost when using this procedure. The important question is then whether a certain metamaterial slab can be treated as a slab of homogeneous material. To answer this question, one can for example calculate the Bloch eigenmodes in an infinitely extended metamaterial and then, in an additional calculation, check which modes are involved when light is reflected by a semi-infinite metamaterial [10,11]. If the calculations show that all but the fundamental Bloch mode are negligible, one can introduce wave parameters, such as the refractive index and wave impedance, for this mode and treat the material as homogeneous. Here, we propose a more straightforward approach, where the properties of a homogenizable metamaterial slab of any thickness are directly linked to the properties of a single layer of the material.
Recently, several retrieval procedures have been introduced to obtain the effective wave parameters, from the reflection and transmission coefficients of a metamaterial slab [10,[12][13][14][15][16][17]. These retrieval procedures rely upon the Fresnel coefficients which are derived for dipolar media. However, for the class of bifacial metamaterials in which the electric quadrupole excitations are present [18], the classical electromagnetic boundary conditions do not hold [19,20]. Consequently, neither the standard Fresnel coefficients nor the retrieval procedures based on these coefficients can be applied to these metamaterials. The development of an adequate theory for the description of highly spatially dispersive metamaterials, such as bifacial metamaterials, would be of great practical importance, e.g. for solar cell applications.
In this work, we re-examine the propagation of light through a slab of spatially dispersive optical metamaterial. We find that both the transmission and reflection by the slab can be surprisingly accurately described in terms of propagating optical plane waves only, which dramatically simplifies the description. We introduce simple analytical expressions that enable one to evaluate the optical characteristics of a thick metamaterial slab by using only the transmission and reflection coefficients of a single layer of the nanoparticles. The singlelayer coefficients are evaluated numerically. These coefficients enable one to calculate the electromagnetic fields inside the material and thereby evaluate the wave parameters, such as the refractive index and wave impedance. Compared to previously reported theoretical approaches to the problem, our approach is easy to use in practice independently of the shapes and material compositions of the nanoparticles and of the propagation direction and polarization of the optical waves as long as the material is homogenizable. Figure 1. A metamaterial slab (a) that is composed by stacking layers of arbitrary scatterers, S, is described as (b) an array of infinitesimally thin sheets. The transmission and reflection coefficients (τ L , τ R , ρ L and ρ R ) for each of these sheets are taken to be equal to those of a single isolated layer of the metamaterial.
Consider a metamaterial slab that is created by stacking two-dimensional arrays of nanoscatterers in a transparent dielectric medium (see figure 1(a)). When illuminated by an optical plane wave, the first layer of nanoscatterers transmits a certain portion of the incident field. A part of this transmitted field is then reflected back by the second layer and the rest is partially absorbed and partially transmitted further to the third layer. Provided that the periodicity within each layer is sufficiently smaller than the illumination wavelength, the diffracted waves will be evanescent and the propagating field between the neighboring layers will consist of two optical plane waves. Furthermore, since the particles in each layer are packed very densely, their collective evanescent field can have a very short decay length in the direction perpendicular to the layer; note that the individual evanescent fields of the particles are still of a long range compared to the particle separation. This effect makes it possible to neglect the nearfields of the layers and describe each layer simply as an infinitesimally thin sheet surrounded by the host dielectric (see figure 1(b)). The total transmission and reflection of the slab is then described in terms of the transmission and reflection coefficients of the individual sheets in a way resembling the description of Fabry-Perot interferometers. Similar division of metamaterial slabs into homogenized layers has been applied in the context of effective electromagnetic parameter retrieval [21,22].
Our approach is as follows. For a single isolated layer of a metamaterial, we first numerically calculate the transmission and reflection coefficients and assign them to an equivalent infinitesimally thin sheet in the middle of the unit cells. These coefficients depend on both the angle of incidence θ and the polarization of the incident field. In addition, for bifacial nanoscatterer arrays, the reflection coefficient changes if the illumination direction is reversed [18,23]. Therefore, for the wave propagating to the right (left) within the metamaterial, we use the reflection coefficient ρ L (ρ R ) to describe the reflection from the left (right) side of each layer. Likewise, the transmission coefficient τ L (τ R ) describes the transmission of a wave incident from the left (right) side of each layer. Optical reciprocity ensures that the transmission coefficients must be the same if normal incidence illumination is considered. When calculated in this way, the parameters ρ L , ρ R , τ L and τ R automatically include the near-and far-field coupling between the scatterers within the layer in question. As has been already mentioned, this coupling makes the extent of the evanescent wave along the layer's normal shorter for denser packing of the scatterers. Let a plane wave with a transverse field amplitude U 0 and a wave vector k in =xk x +ẑk z be incident on a metamaterial slab that consists of N layers of thickness and has its surface normal along the z-axis (see figure 2). The wave is assumed to have either transverse electric (TE) or transverse magnetic (TM) polarization. Treating each nanoscatterer layer as an infinitesimally thin sheet, we consider the counter-propagating waves between the sheets to have wave vectors k ± =xk x ±ẑk z , because the material between the sheets is considered to be the same as outside the slab. Assuming that the polarization state is conserved, the transverse fields U j and U j after each sheet j in figure 2 must satisfy the following equations: For chiral metamaterials, both polarization states must be considered simultaneously. In this case, U j and U j would be two-element vectors, while f L , f R , g L and g R would be 2 × 2 matrices.
Using (1) and (2), we derive separate relations for the forward and backward propagating fields: where With the help of (1) and (3) and the fact that U N = 0, we obtain the transmission coefficient of the slab to be Here we have introduced the G-polynomial that is calculated as Similarly, using (1), (2) and (4) we derive the reflection coefficient Equations (5) and (9) enable direct calculation of the transmission and reflection coefficients of an arbitrarily thick metamaterial in terms of the transmission and reflection coefficients of an isolated monolayer of the metamaterial. For N = 1, equations (5) and (9) correctly yield t = τ L exp(ik z ) and r = ρ L exp(ik z ) and for N = 2 the well-known results for a Fabry-Perot etalon are obtained. In order to demonstrate the applicability of our theory, we compare it with rigorous numerical calculations. This is done by selecting some non-trivial nanoscatterers, such as nanoshells (figure 3(a)), nanorings (figure 4(a)) and nanodimers ( figure 5(a)). These nanoscatterers are considered to compose stacks of two-dimensional periodic arrays that are embedded in a dielectric host medium of refractive index 1.5. The necessary transmission and reflection coefficients for a single array are calculated using the computer software COMSOL Multiphysics. The obtained coefficients are then used in (5) and (9) to acquire the transmission and reflection coefficients for several layers. These coefficients are compared with the results of direct numerical calculations of the whole stack with COMSOL. For these calculations we choose TM-polarized illumination with θ = 45 • and a slab consisting of N = 5 layers. This choice is general enough for demonstrating the applicability of the model.
We first consider a metamaterial with an isotropic unit cell containing a silver nanoshell as depicted in figure 3(a). The nanoshells have an outer radius R = 35 nm and a thickness h = 7 nm. They form a cubic lattice with period = 130 nm. The optical properties of such nanoshells are well investigated [24][25][26] and similar structures can be relatively readily fabricated [27,28]. For the calculations, the optical characteristics of silver were taken from [29]. The calculated transmission and reflection spectra in the wavelength range from 500 to 1000 nm for the nanoshell slab are shown in figure 3(b). The spectra obtained by using (5) and (9) are in a very good agreement with the direct numerical calculations, indicating that the plane-wave description of the light propagation is appropriate and that the evanescentwave coupling between the layers is indeed weak. A more detailed description of the influence of this coupling on the transmission and reflection spectra is presented later on in the paper.
While direct numerical calculations are always limited by the present computational resources, the introduced theory allows us to calculate the response of an arbitrarily thick slab. In figure 3(c) the reflection spectrum for a slab of 1000 nanoshell layers is shown by the red curve. For this particular case of nanoshells, the spectrum of 1000 layers is close to the spectrum of   (5) and (9). (c) Reflection spectra for increasing number of layers as obtained from (5) and (9). the five-layer slab and it is already indistinguishable from the spectrum of an infinitely thick metamaterial.
Next, we introduce an anisotropic (uniaxial) unit cell containing a silver nanoring as depicted in figure 4(a). An interesting application of such structures as optical security marks is proposed in [30]. The ring has an outer radius R 1 = 20 nm, inner radius R 2 = 10 nm and thickness h = 10 nm. The rings form a cubic lattice with period = 50 nm, such that each layer is aligned with the xy-plane. Note that within each layer the interparticle separation distance is only 10 nm. The calculated transmission and reflection spectra for the nanoring slab are shown in figure 4(b). The theory yields excellent agreement with direct numerical calculations also for these nanoscatterers. The reflection for 1000 layers, depicted by the red curve in figure 4(c),  shows that the bulk metamaterial behaves quite differently from a single layer due to the interlayer interaction. One can notice that if the number of layers is large, the metamaterial acts as a spectrally selective broad-band reflector with a nearly flat-top spectral profile. As a final example, we consider a bifacial metamaterial slab that exhibits strong spatial dispersion. Such metamaterials have not been studied much in terms of their reflection and transmission characteristics. The unit cells of the material contain asymmetric silver nanodimers (see figure 5(a)). These nanodimers have been shown to exhibit complete suppression of the electric dipole excitation in a narrow wavelength range when illuminated from the smaller disc side [18]. However, in a metamaterial, there will be two counter-propagating waves and the electric dipole moment cannot be suppressed for both of them simultaneously [23]. The nanodimer geometry is described by the radii R 1 = 15 nm and R 2 = 20 nm of the discs and dimensions h = s = 10 nm defined in figure 5(a). A cubic lattice with period = 50 nm is now composed of the nanodimers such that the smaller discs are on the left-hand side. As a consequence of the asymmetry of the unit cell, we calculate the single layer response to  figure 3(b). The numerically and analytically calculated reflections from the side of the larger discs are shown by the additional red dotted line and red triangles, respectively. (c) Reflection spectra for 1000 layers as obtained from (5) and (9).
illumination from both sides in order to obtain the reflection coefficients ρ L and ρ R . The invariance of the nanodimers with respect to rotation around the symmetry axis z ensures that τ L = τ R . Using (5) and (9) we calculate the transmission and reflection spectra for a nanodimer slab illuminated from the two sides and compare them with the numerical results. Figure 5(b) shows that while the theory very accurately resolves all spectral features, there is a slight deviation of the analytically obtained values from the exact numerical values for the wavelengths around 600 nm. This deviation obviously originates from the evanescent-wave coupling between the adjacent layers. However, considering that the gap size between the discs in the adjacent layers is only 20 nm, the agreement is still remarkably good. We obtained a similar good agreement between our theory and the numerical calculations also for angles of incidence of 0 • , 30 • and 60 • as well as for the TE polarization. The reflection spectra of 1000 layers of nanodimers are depicted in figure 5(c). When illuminated from the small disc side, the reflection coefficient significantly decreases at around the electric dipole suppression wavelength of 618 nm.
Using (1)-(4), one can retrieve the electric and magnetic fields at any point inside the metamaterial and use them to directly extract the effective wave parameters, such as the refractive index and wave impedance. As an example, consider a non-chiral centrosymmetric material, for which f L = f R = f and g L = g R = g. Propagation of a plane wave over a single unit cell in the homogenized material must satisfy where γ z is the z component of the effective propagation constant and is the unit-cell size in the z direction. Equation (3) then leads to the following expression: where m is an integer. The effective refractive index is related to the wave vector k in the host medium through which follows from the phase matching condition k x = γ x ; k 0 is the wave number in vacuum. The effective wave impedance can be obtained by considering the spatially averaged electric and magnetic fields between the sheets introduced in figure 1. Since in the host medium the fields are right-handed, the total electric field component that is transverse to z is U j exp(ik z z) + U j exp(−ik z z), whereas the magnetic field component that is transverse to z is [U j exp(ik z z) − U j exp(−ik z z)][k z /k] p /Z , with p = ±1 denoting the TE and TM polarizations, respectively. Here Z denotes the wave impedance in the host medium. We can now define that describes the ratio between the transverse components of the averaged electric and magnetic fields; the angle brackets denote averaging over the unit cell. Taking into account the fact that the propagation angle of the effective wave is determined by γ z /γ , one can find the effective wave impedance Z eff = Z ⊥ (γ z /γ ) p . Performing the averaging in (13), we obtain where U j was expressed in terms of U j using (1) and (10). We recall that γ = n eff k 0 . One can also obtain the corresponding relative electric permittivity and magnetic permeability as  where Z 0 is the wave impedance in vacuum. It can be verified that the expressions for n eff , Z eff , ε eff and µ eff are in full agreement with the commonly used expressions introduced in [15]. The derivations above can be repeated also for more complex materials, with f L = f R and g L = g R . The parameters n eff , Z eff , ε eff and µ eff calculated for the nanoring material of figure 4 are shown in figure 6, for θ = 45 • and TE polarization. This example is of interest in view of the possibility to tune n eff and Z eff , because the rings are somewhat similar to traditional split-ring resonators. A strong electric-dipole resonance at around λ 0 = 870 nm is observed in the spectra of n eff and ε eff . The modification of µ eff in this spectral range is not large and can be interpreted as a result of the finite periodicity (see e.g. [31]). At wavelengths shorter than λ 0 = 800 nm, the behavior of ε eff resembles that of a Drude metal. A thick nanoring material could therefore have a high reflectivity in this region, which is supported by the values of the wave impedance. It is interesting, however, that at λ 0 ≈ 590 nm the material is characterized by n eff ≈ 1 and Z eff ≈ Z 0 , leading to an efficient suppression of both optical reflection and refraction at an air-metamaterial interface.
The accuracy of the presented theory depends on the extent of the evanescent waves produced by the nanoparticle layers. Qualitatively, for the theory to be exact, the evanescent waves associated with the cut-off diffraction orders must have a decay length that is much shorter than the spacing between the particles in two adjacent layers. For a two-dimensional square array of period x = y , the longitudinal wave vector of the first such order k z1 = [k 2 − (2π/ x − k x ) 2 ] 1/2 is imaginary, with k and k x being the magnitudes of the total and transverse wave numbers of the incident light. The next order would be k z2 = [k 2 − (2π/ x + k x ) 2 ] 1/2 . We define the decay length δ of the evanescent field to be the distance for which the field amplitude   (5) and (9).
has decayed by a factor of exp(−1). The gap d between the particles in the z-direction must then be much larger than Considering the nanoshells with x = 130 nm, λ 0 = 500 nm and k x = k/ √ 2, we obtain a decay length δ ≈ 34 nm that is smaller than the 60 nm gap between the adjacent shells. For the nanodimers, with x = 50 nm, λ 0 = 500 nm and k x = k/ √ 2, the decay length is δ ≈ 9 nm. This value is smaller than the gap d = 20 nm between the nanodimers, which supports the success of our analytical calculations. In fact, if for subwavelength-sized unit cells we have k 2π/ x , (17) yields δ ≈ x /(2π ). In this case the criterion for neglecting the interlayer evanescent-wave coupling becomes d x /(2π). Then, as a practical criterion for when our theory can be applied, we require that d > x /2.
In order to verify the above predictions on the influence of the interlayer evanescent-wave coupling, we numerically calculate the transmission through the nanodimer slab, while varying the transverse and longitudinal periods separately. In figure 7(a) the transmission coefficient is plotted for an increasing longitudinal period z . The transverse period is fixed to x = 100 nm in order to have the evanescent-wave coupling significant enough to cause a discrepancy at small z between the theory and the numerical results. We notice that when d exceeds x /2 ( z exceeds 80 nm), this discrepancy disappears. The transmission coefficient for an increasing transverse period x = y is shown in figure 7(b). In this case the longitudinal period is fixed to z = 70 nm, such that the evanescent-wave coupling is negligible at x = 50 nm. As the transverse period is increased, the discrepancy between the theory and the numerical results starts to appear due to an increase in the evanescent-wave coupling between the layers. We notice however that as x is increased, the array also gets sparse, which reduces the influence of the nanoscatterers on the propagating wave. This effect counterbalances the growing decay length of the evanescent waves produced by the nanoparticle layers, such that the discrepancy between the theory and the numerical results remains small.
For a metamaterial slab, in which the evanescent wave coupling is negligible, one can introduce effective material parameters. On the other hand, if the evanescent-wave coupling exists, these parameters depend on the slab thickness and are thereby senseless [32]. By using (5) and (9) to compare the transmission and reflection coefficients of a single nanoscatterer layer with those of two layers, one can directly assess whether the metamaterial is homogenizable and, consequently, whether the introduction of material parameters is justified.
In summary, we have introduced a simple analytical theory for the description of light interaction with optical metamaterials. Recognizing the subwavelength size of the metamaterial's unit cells, we found that the evanescent-wave coupling between adjacent monolayers of the metamaterial does not significantly influence the light propagation in the material. As a consequence, arbitrarily thick metamaterial slabs can be accurately described in terms of the plane-wave transmission characteristics obtained for a single isolated monolayer, e.g. numerically. Furthermore, we have shown that one can evaluate the fields at any point inside the material and, consequently, obtain the effective wave parameters. The presented examples of rigorous numerical calculations demonstrate the wide applicability of this remarkably simple analytical model.
In contrast to existing theoretical approaches, our one also correctly describes threedimensional arrays of bifacial nanoscatterers, which is of practical importance for a large variety of metamaterials, such as those with asymmetric unit cells. For homogenizable metamaterial slabs, our method enables rapid one-layer-based extraction of the transmission and reflection coefficients. Furthermore, propagation of an optical beam through such a metamaterial can be described by using the angular spectrum representation with our model applied to each planewave component.
The introduced theory is not limited to optical metamaterials, but can also be applied to study wave propagation in other artificial media, such as radio-frequency and terahertz metamaterials, and even phononic metamaterials. We believe that the presented theory has the necessary simplicity and accuracy to accelerate the development of optical metamaterials tailored for real applications.