Ultrathin metamaterial layers are modeled by a homogeneous bi-anisotropic film to represent various kinds of broken symmetries in photonic nanostructures, and specifically in optical metamaterials and metasurfaces. Two algorithms were developed to obtain the electromagnetic (EM) wave response from a metasurface (direct solver) or the metasurface parameters from the EM wave response (inverse solver) for a bi-anisotropic, subwavelength-thick nanostructured film. The algorithm is applied to two different metasurfaces to retrieve their effective homogeneous bi-anisotropic parameters. The effective layer of the same physical thickness is shown to produce the same response to plane wave excitation as the original metasurface.
© 2013 Optical Society of America
Metamaterials have been used for manipulating light in a controllable manner, and they were able to achieve optical properties not existing in nature including negative index of refraction [1, 2], optical magnetism , invisibility cloaking  and superlensing [5, 6]. Optically thin metamaterial layers called optical metasurfaces have been introduced later, and they have the ability to control light through engineering of the wavefront of the incident light. Metasurfaces have garnered interest due to the fact that they avoid the significant losses associated with most bulk metamaterials, and are suitable for on-chip application, and also their elements are more readily assembled . Metasurfaces have been used to implement important applications such as light bending [8, 9], flat lenses , circular polarizers , half-wave plates [12, 13], and quarter wave plates [14, 15]. Those metasurfaces provide their intended functionality by changing the phase and/or polarization of light transmitted through the layer. In order to develop metasurfaces and fully analyze their functionalities, it is important to have an accurate and efficient model to describe the unit cell of the surface nanostructure. In this work, we have developed a model for metasurface layers with a thin, homogeneous, equivalent film. Using this framework, metasurface designers can then obtain insight on how best to use the unit-cell structures.
Most of metasurface designs depend on symmetry breaking in the nanostructure, such as rotational symmetry, mirror symmetry or directional symmetry. A bi-anisotropic model would be quite general and useful option for the homogenization of metasurface designs. Throughout the paper, the part of the developed model responsible for each kind of asymmetry will be explained. The goal is to obtain a homogenous bi-anisotropic film that will generate the same values of the complex coefficients for reflection and transmission as those obtained by the real metamaterial structure. Homogenization using reflection and transmission coefficients has been used for permittivity and permeability retrieval , but non-physical dispersion relations may occur due to limitation of the model, and lack of representation of bi-anisotropy, chirality and spatial dispersion . Bi-anisotropy has been used in  to account for directional asymmetry, and in this work a general bi-anisotropic tensor is used to account for bi-anisotropy, chirality and linear spatial dispersion. The model is only limited by higher order spatial dispersion terms, however, their effect will significantly decrease in thin fillms and the linear spatial dispersion term will dominate. Detailed discussions on the physical meaning of metamaterial constitutive parameters can be found in [17, 19]. The work is done by first solving for the transmission and reflection coefficients of a bi-anisotropic layer developing a direct solver. Then, it is solved for a bi-anisotropic film that accurately represents our structure, thus developing an inverse solver. After explaining the details of the process, the algorithm is implemented to homogenize two specific structures that are commonly used in metasurface applications.
2. Direct solver
The first step in the homogenization process is to compute the complex reflection and transmission coefficients for a thin, bi-anisotropic film of known parameters. Table 1 clarifies the complex reflection and transmission coefficients (4-port S-parameters).
As shown in Table 1 and further illustrated in Fig. 1, the subscripts of S-parameters represent the output and input sides respectively (1 – superstrate, 2 – substrate), while their superscripts represent the polarization of the output and input waves respectively. There also exist the other set of eight coefficients for the waves incident at the substrate side (side 2).
The co-polarized reflection relation is described by using the complex co-polarized transmission coefficient , where . In contrast, the cross-polarized reflection relation is described by using the complex co-polarized reflection coefficient , where for . In the both cases (of co- or cross-polarized) reflections , for ‘ + ’, or , for ‘-’. The co-polarized transmission relation is described by , where . The cross-polarized transmission is defined by , where for . In both cases (either co- or cross-polarized) reflections , for ‘ + ’, or , for ‘-’.
The electromagnetic (EM) waves propagating inside the bi-anisotropic structure satisfy the following material equation for the field components of a normally incident plane wave:20], and we cannot decouple the transverse-electric (TE) and the transverse-magnetic (TM) waves inside the film.
To solve for the eigenmodes of the system, we start from Maxwell’s equations:Equation (2) can be formally rewritten as
The eigendecomposition of the transfer matrix yields
The eigenvectors , ,, represent four possible wavefronts propagating through the bi-anisotropic medium, and the eigenvalues are the wavenumbers corresponding to each eigenmode. Typically, two of the eigenvalues are positive and two are negative corresponding to forward and backward propagation, respectively. Thus, in which combines (4) and (5), the term decomposes as a superposition of the four eigenvectors. Then the matrix applies a propagation term to each eigenvector. Finally, the matrix sums all four propagated eigenvector components of .
The reflection and transmission coefficients are obtained as we apply an x-polarized and a y-polarized wave. First, we apply an x-polarized input wave at the front-side (superstrate- side) with the normalized values and , where and are respectively the intrinsic impedances of the superstrate and substrate. Then, we have the input-side vector due to the addition of incident and reflected wave in the following form , , and the output vectors being , , where , and and are the complex reflection and transmission coefficients.
Therefore, after partitioning the transfer matrix ( is the thickness of the BA layer) and writing
In a similar way, for y-polarized input , we can obtain the complex reflection and transmission coefficients for y-polarized superstrate-side illumination:
For back illumination (substrate illumination), using same route, we get , , with the output vectors being , .Using (6), we get:
Similarly for y-polarized input, we get:
And hence, all the complex reflection and transmission coefficients described in Table 1 are obtained. One major advantage of this direct solver is that it depends on simple matrix operations, which are reversible. This makes the development of the inverse solver straightforward as described in the following section.
3. Inverse problem
The inverse problem uses the complex reflection and transmission coefficients to obtain the corresponding material dyadics described in Eq. (1). This is accomplished in two steps. First the linear operator is retrieved; then it is used to obtain all the material constants. For x- and y-polarized inputs used at the front ‘ + ’ and back ‘‒‘ side illumination, there are four sets of fields , that can be used to form the equation:
Moreover, as we recall from (3) that
This concludes the retrieval of the BA parameters from a given set of complex reflection and transmission coefficients obtained upon four distinct illumination states.
It is important to note that the eigendecomposition of , which is where is the diagonal matrix carrying the eigenvalues of , could suffer from phase ambiguity if one of the phase terms which are real parts of , , or is above or below , This should not be the case for a metasurface of sub-wavelength thickness. Some techniques have been developed to overcome phase ambiguity for the retrieval of bulk media such as the one in  which performs the retrieval over a spectral range where there is no ambiguity at the largest wavelength. The technique then utilizes phase unwrapping along the remaining spectral range to remove the ambuguity. This can be implemented to our technique by applying phase unwrapping to the elements of the diagonal matrix , but we don’t need it in our work with metasurfaces.
It can be confirmed that the formula in Eq. (20) converges to the analytical formulas obtained for special cases of material properties. The simplest of them is the isotropic medium with properties,,,surrounded by free space of impedance . For this case, due to symmetry, we have the same co-polarized complex reflection and co-polarized transmission for all illumination cases and zero cross-polarization coefficients. Therefore, we have:
The eigendecomposition of the matrix in (22) is required to apply the formula in (20). Solving for the eigenvalues will yield the following equation using the symbol for eigenvalues:
It’s clear that these eigenmodes correspond to forward and backward x- and y-polarized plane waves. Applying the formula in (20):16]. The next step is to apply the inverse problem described by (20) to some metasurface structure and explore how different kind of broken symmetries are represented in the effective material model.
4. Application of algorithm to metasurfaces
In this section, the algorithm is implemented to homogenize two different metasurfaces. One metasurface is an array of V-shaped, gold antennas with a thickness of 30 nm fabricated on top of a silicon substrate. Figure 2(a) a unit-cell of structure with dimensions of 200 nm x 200 nm, and a V-shape angle of 60° between the two arms. Each arm has a length of 160 nm, and a width of 40 nm. The second structure, shown schematically in Fig. 2(b), has a unit-cell size of 300 nm x 300 nm and is composed of 2 gold rods, each of a 250-nm length, a 40-nm width, a 30-nm thickness, a vertical separation of 80 nm, and an orientation angle of 45°. These rods are embedded in a 200-nm-thick polymer layer. Both structures will have a broken rotational symmetry (anisotropy).
The diagonal terms of and cause coupling between x- and y-polarized waves as demonstrated in the studies done to the bi-isotropic case , while their off-diagonal terms affect the relation between the electric and magnetic fields while keeping x- or y-polarization, but causing the waves to experience different wave impedances for propagation in + z or –z direction. The structure of Fig. 2(a) also has a broken directional symmetry due to the difference between superstrate and substrate plasmonic resonances, but the structure keeps its mirror symmetry with respect to x-axis causing negligible coupling between x- and y-polarized waves. Therefore, the model of this structure should contain only the off-diagonal terms of and and the diagonal terms should be negligible. However, the structure in Fig. 2(b) lacks mirror symmetry with respect to x- or y-axis but has directional symmetry with respect to z-direction. Therefore its model should contain only the diagonal terms of and . The inverse solvers were applied to transmission and reflection coefficients of both structures (obtained using FEM) and the results were as expected. Figure 3 shows the results of implementing the homogenization algorithm with the V-antenna structure.
The V-antenna structure is modeled as a 30-nm-thick, homogeneous slab. First, the FEM solver is used to obtain the complex reflection and transmission coefficients for a spectral domain of 1µm – 3µm. Then the inverse solver is applied to obtain the effective parameters of the slab. The retrieved results of the effective slab parameters are shown in Fig. 3(a). This design has been used to achieve large phase shifts [8, 9] from operation near the plasmonic resonance wavelength, and indeed the effect of the resonance is clear in the retrieved effective parameters. Both the FEM simulation and the effective bi-anisotropic model reproduce the same reflection and transmission coefficients as shown in Fig. 3(b).
In this structure, we obtain only co-polarized reflection and transmission coefficients, with negligible values for cross-polarization coefficients. The directional asymmetry is noticed in the difference between reflection coefficients for the two sides of illumination. Still the transmission coefficients are symmetric (i.e. and , where is the refractive index of the backward substrate ), and this would result in having and as mathematically proved in , so we needed only to show and in Fig. 3(a).
Now the algorithm is applied to the two rod structure which is modeled as a homogeneous bi-anisotropic slab with a thickness of 200 nm. The FEM solver is used to obtain the complex reflection and transmission coefficients for the spectral domain of 2µm – 3µm. The retrieval algorithm is then applied to these data, and as in the previous case, the complex reflection and transmission coefficients obtained from the effective model match with these obtained from FEM simulation. The retrieval results are shown below in Fig. 4. This structure has a mirror asymmetry or parity asymmetry and this results in the existence of the diagonal elements of the tensors and .This structure is exactly the same from both sides, and this directional symmetry causes the off-diagonal terms of and to be zero.
In this paper, a new approach has been presented for the homogenization of optical metamaterials. An algorithm was developed which included a direct and an inverse solver based on an eigenwave analysis and a transfer matrix approach. This method of modeling and characterizing a metamaterial is useful for the design and use of metasurfaces. In addition, it provides insight into how these metamaterial layers affect an incident light beam. The used model is the most general bi-anisotropic model for normal incident illumination. This model can be extended to include effects on oblique incident waves. .
This work was partially supported by AFOSR grant FA9550-10-1-0264, ARO grant 57981-PH (W911NF-11-1-0359), and NSF grant DMR-1120923.
References and links
1. V. M. Shalaev, “Optical negative-index metamaterials,” Nat. Photonics 1(1), 41–48 (2007). [CrossRef]
3. U. K. Chettiar, S. Xiao, A. V. Kildishev, W. Cai, H. K. Yuan, V. P. Drachev, and V. M. Shalaev, “Optical metamagnetism and negative-index metamaterials,” MRS Bull. 33(10), 921–926 (2008). [CrossRef]
4. D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science 314(5801), 977–980 (2006). [CrossRef] [PubMed]
8. N. Yu, P. Genevet, M. A. Kats, F. Aieta, J. P. Tetienne, F. Capasso, and Z. Gaburro, “Light propagation with phase discontinuities: generalized laws of reflection and refraction,” Science 334(6054), 333–337 (2011). [CrossRef] [PubMed]
10. F. Aieta, P. Genevet, M. A. Kats, N. Yu, R. Blanchard, Z. Gaburro, and F. Capasso, “Aberration-free ultrathin flat lenses and axicons at telecom wavelengths based on plasmonic metasurfaces,” Nano Lett. 12(9), 4932–4936 (2012). [CrossRef] [PubMed]
12. J. Hao, Y. Yuan, L. Ran, T. Jiang, J. A. Kong, C. T. Chan, and L. Zhou, “Manipulating electromagnetic wave polarizations by anisotropic metamaterials,” Phys. Rev. Lett. 99(6), 063908 (2007). [CrossRef] [PubMed]
13. Z. H. Zhu, C. C. Guo, K. Liu, W. M. Ye, X. D. Yuan, B. Yang, and T. Ma, “Metallic nanofilm half-wave plate based on magnetic plasmon resonance,” Opt. Lett. 37(4), 698–700 (2012). [CrossRef] [PubMed]
15. Y. Zhao and A. Alù, “Manipulating light polarization with ultrathin plasmonic metasurfaces,” Phys. Rev. B 84(20), 205428 (2011). [CrossRef]
16. D. Smith, S. Schultz, P. Markoš, and C. Soukoulis, “Determination of effective permittivity and permeability of metamaterials from reflection and transmission coefficients,” Phys. Rev. B 65(19), 195104 (2002). [CrossRef]
17. A. Alù, “First-principles homogenization theory for periodic metamaterials,” Phys. Rev. B 84(7), 075153 (2011). [CrossRef]
18. A. V. Kildishev, J. D. Borneman, X. Ni, V. M. Shalaev, and V. P. Drachev, “Bianisotropic effective parameters of optical metamagnetics and negative-index materials,” Proc. IEEE 99(10), 1691–1700 (2011). [CrossRef]
19. A. Alù, “Restoring the physical meaning of metamaterial constitutive parameters,” Phys. Rev. B 83(8), 081102 (2011). [CrossRef]
20. I. V. Lindell and A. J. Viitanen, “Eigenwaves in the general uniaxial bianisotropic medium with symmetric parameter dyadics,” Report No. 148 Helsinki Univ. of Technology, Espoo (Finland). Electromagnetics Lab. 1 (1993).
21. O. Luukkonen, S. I. Maslovski, and S. A. Tretyakov, “A stepwise Nicolson–Ross–Weir-based material parameter extraction method,” IEEE Antennas Wirel. Propag. Lett. 10, 1295–1298 (2011). [CrossRef]
22. D.-H. Kwon, D. H. Werner, A. V. Kildishev, and V. M. Shalaev, “Material parameter retrieval procedure for general bi-isotropic metamaterials and its application to optical chiral negative-index metamaterial design,” Opt. Express 16(16), 11822–11829 (2008). [CrossRef] [PubMed]