Homogenization of bi-anisotropic metasurfaces

Ultrathin metamaterial layers are modeled by a homogeneous bianisotropic 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 bianisotropic 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 OCIS codes: (160.3918) Metamaterials; (310.0310) Thin Films; (100.3190) Inverse Problems. References and links 1. V. M. Shalaev, “Optical negative-index metamaterials,” Nat. Photonics 1(1), 41–48 (2007). 2. R. A. Shelby, D. R. Smith, and S. Schultz, “Experimental verification of a negative index of refraction,” Science 292(5514), 77–79 (2001). 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). 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). 5. J. B. Pendry, “Negative refraction makes a perfect lens,” Phys. Rev. Lett. 85(18), 3966–3969 (2000). 6. N. Fang, H. Lee, C. Sun, and X. Zhang, “Sub-diffraction-limited optical imaging with a silver superlens,” Science 308(5721), 534–537 (2005). 7. N. Engheta, “Antenna-guided light,” Science 334(6054), 317–318 (2011). 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). 9. X. Ni, N. K. Emani, A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, “Broadband light bending with plasmonic nanoantennas,” Science 335(6067), 427–427 (2012). 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). 11. Y. Zhao, M. A. Belkin, and A. Alù, “Twisted optical metamaterials for planarized ultrathin broadband circular polarizers,” Nat. Commun. 3, 870 (2012). 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). 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). 14. A. Drezet, C. Genet, and T. W. Ebbesen, “Miniature plasmonic wave plates,” Phys. Rev. Lett. 101(4), 043902 (2008). 15. Y. Zhao and A. Alù, “Manipulating light polarization with ultrathin plasmonic metasurfaces,” Phys. Rev. B 84(20), 205428 (2011). 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). 17. A. Alù, “First-principles homogenization theory for periodic metamaterials,” Phys. Rev. B 84(7), 075153 (2011). 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). 19. A. Alù, “Restoring the physical meaning of metamaterial constitutive parameters,” Phys. Rev. B 83(8), 081102 (2011). #185480 $15.00 USD Received 4 Mar 2013; revised 27 Jun 2013; accepted 13 Aug 2013; published 11 Sep 2013 (C) 2013 OSA 23 September 2013 | Vol. 21, No. 19 | DOI:10.1364/OE.21.021941 | OPTICS EXPRESS 21941 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). 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).


Introduction
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 [3], invisibility cloaking [4] 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 [7].Metasurfaces have been used to implement important applications such as light bending [8,9], flat lenses [10], circular polarizers [11], 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 [16], but non-physical dispersion relations may occur due to limitation of the model, and lack of representation of bianisotropy, chirality and spatial dispersion [17].Bi-anisotropy has been used in [18] 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.

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

S
As shown in Table 1 and further illustrated in Fig. 1, the subscripts of S-parameters represent the output and input sides respectively (1superstrate, 2substrate), 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 qq q q jj j j S E E   using the complex copolarized transmission coefficient qq jj S , where   x, y q  . In contrast, the cross-polarized reflection relation is described by pq q p jj j j S E E   using the complex co-polarized reflection coefficient pq jj S , where   x, y p  for   y, x q  .In the both cases (of co-or cross-polarized) reflections 1 j  , for ' + ', or 2 j  , for '-'.The co-polarized transmission relation is described by qq q q ij j i

S E E  
, where   x, y q  . The cross-polarized transmission is defined by , where

 
x, y p  for   y, x q  . In both cases (either co-or cross-polarized) reflections 1 j  , for ' + ', or 2 j  , 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: .The parameters ε, are the relative permittivity and relative permeability tensors, and ξ, are the bi-anisotropic tensors.The free space constants of permittivity 0  , permeability 0  and speed of light c are used to normalize the model.We use tensors to represent the rotational asymmetry (anisotropy).The diagonal terms of ξ,ζ model the mirror asymmetry effect (chirality) and their off-diagonal terms are responsible for breaking the directional symmetry of the propagating wave.The eigenmodes excited inside the bi-anisotropic medium are in general elliptically polarized [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:  2) can be formally rewritten as where

AN
, where , and 01 10 ; it has a straightforward solution of: which is used to directly obtain the field components at any location inside the biansiotropic media using the components at the origin where V 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 4) and ( 5), the term 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 (superstrateside) with the normalized values x,inc 1 E  and y,inc , where 1 z and 2 z 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 I S , and the output vectors being T n , and In a similar way, for y-polarized input

S
coefficients for y-polarized superstrate-side illumination: For back illumination (substrate illumination), using same route, we get .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.

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 T 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: where all four illumination states are grouped together as:
From ( 17) and the connection between T and A given by exp l 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 k 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 [21] 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 P , 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 0 ε= 0 , ξ=0 , = 0 surrounded by free space of impedance 0 0 0 z   .For this case, due to symmetry, we have the same co-polarized complex reflection r and co-polarized transmission t for all illumination cases and zero cross- polarization coefficients.Therefore, we have: where 10 01

I
, and T is given by:  (22) 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 p for eigenvalues:   or, and solving for the corresponding eigenvectors, we get: It's clear that these eigenmodes correspond to forward and backward xand y-polarized plane waves.Applying the formula in (20): which leads to: where k and z are given by ( 25) and ( 27).The previous results are identical to retrieval formulas in [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.

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 [22], 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 ypolarized 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 .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. 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.

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

Fig. 1 .
Fig. 1.Demonstration of co-polarized and cross-polarized reflection and transmission coefficients for an x-polarized input plane wave from the front size.The incident wave is indicated in blue, and the reflected and transmitted waves are in red.The direction of (E) and (H) fields are shown by arrows.
P is the diagonal matrix carrying the eigenvalues of T , could suffer from phase ambiguity if one of the phase terms which are real parts of a kl, b kl, c l k or d l

Fig. 2 .
Fig. 2. The unit cell of the two nanostructures used as examples for the homogenization algorithm.A top view is presented for each to describe the x and y directions.

Fig. 3 .
Fig. 3. (a) Effective slab parameters of the V-antenna structure.(b) Spectral phasor of reflection and transmission coefficients using FEM and effective Bi-anisotropic layer model.For all of them, the lower end point corresponds to 1 µm and the upper end point corresponds to 3 µm.

Fig. 4 .
Fig. 4. Retrieval results of the two rod structure.

Table 1 . Complex reflection and transmission coefficients
