1. 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 bi-anisotropy, 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.

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

Table 1. Complex reflection and transmission coefficients

View Table

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

 

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.

Download Full Size | PDF

The co-polarized reflection relation is described by $S_{jj}^{qq}{E}_{j\pm }^q=E_{j\mp }^q$ using the complex co-polarized transmission coefficient $S_{jj}^{qq}$, where $q=\left\{\text{x},\text{y}\right\}$. In contrast, the cross-polarized reflection relation is described by $S_{jj}^{pq}{E}_{j\pm }^q=E_{j\mp }^p$ using the complex co-polarized reflection coefficient $S_{jj}^{pq}$, where $p=\left\{\text{x},\text{y}\right\}$ for $q=\left\{\text{y},\text{x}\right\}$. In the both cases (of co- or cross-polarized) reflections $j=1$, for ‘ + ’, or $j=2$, for ‘-’. The co-polarized transmission relation is described by $S_{ij}^{qq}{E}_{j\pm }^q=E_{i\mp }^q$, where $q=\left\{\text{x},\text{y}\right\}$. The cross-polarized transmission is defined by $S_{jj}^{pq}{E}_{j\pm }^q=E_{j\mp }^p$, where $p=\left\{\text{x},\text{y}\right\}$ for $q=\left\{\text{y},\text{x}\right\}$. In both cases (either co- or cross-polarized) reflections $j=1$, for ‘ + ’, or $j=2$, 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:

To solve for the eigenmodes of the system, we start from Maxwell’s equations:

The eigendecomposition of the transfer matrix $T\left(z\right)=e^{\iota Az}$ yields

The eigenvectors $V_{\text{a}}$, $V_{\text{b}}$,$V_{\text{c}}$,$V_{\text{d}}$ 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 $V\left(z\right)=UP{U}^{-1}{V|}_{z=0}$ which combines (4) and (5), the term $U^{-1}$ decomposes ${V|}_{z=0}$ as a superposition of the four eigenvectors. Then the matrix $P$ applies a propagation term to each eigenvector. Finally, the matrix $U$ sums all four propagated eigenvector components of $V\left(z\right)$.

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 $E_{\text{x,inc}}=1$ and $H_{\text{y,inc}}=1/z_1$, where $z_1$and $z_2$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 $E_1=I_{\text{x}}+S_{11}^{\text{x}}$, $H_1=z_1^{-1}n\left(I_{\text{x}}-S_{11}^{\text{x}}\right)$, and the output vectors being $E_2=S_{21}^{\text{x}}$, $H_2=z_2^{-1}n{S}_{21}^{\text{x}}$, where $I_{\text{x}}={\left(\begin{array}{cc}1& 0\end{array}\right)}^{\text{T}}$, and $S_{11}^{\text{x}}={\left(\begin{array}{cc}{S}_{11}^{\text{xx}}& S_{11}^{\text{yx}}\end{array}\right)}^{\text{T}}$ and $S_{21}^{\text{x}}={\left(\begin{array}{cc}{S}_{21}^{\text{xx}}& S_{21}^{\text{yx}}\end{array}\right)}^{\text{T}}$ are the complex reflection and transmission coefficients.

Therefore, after partitioning the transfer matrix $T\left(l\right)=\exp \iota Al=\left(\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right)$ ($l$ is the thickness of the BA layer) and writing

In a similar way, for y-polarized input $I_{\text{y}}={\left(\begin{array}{cc}0& 1\end{array}\right)}^{\text{T}}$, we can obtain the complex reflection $S_{11}^{\text{y}}={\left(\begin{array}{cc}{S}_{11}^{\text{xy}}& S_{11}^{\text{yy}}\end{array}\right)}^{\text{T}}$ and transmission $S_{21}^{\text{y}}={\left(\begin{array}{cc}{S}_{21}^{\text{xy}}& S_{21}^{\text{yy}}\end{array}\right)}^{\text{T}}$ coefficients for y-polarized superstrate-side illumination:

For back illumination (substrate illumination), using same route, we get $E_2=I_{\text{x}}+S_{22}^{\text{x}}$, $H_2=z_2^{-1}n\left(S_{22}^{\text{x}}-I_{\text{x}}\right)$, with the output vectors being $E_1=S_{12}^{\text{x}}$, $H_1=-z_1^{-1}n{S}_{12}^{\text{x}}$.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 $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 $V_j=\left(\begin{array}{cc}{V}_j^{+}& V_j^{-}\end{array}\right)$, $j=\left\{1,2\right\}$ that can be used to form the equation:

Moreover, as we recall from (3) that $A=\omega N\left(\begin{array}{cc}{\epsilon }_0\text{ε}& c^{-1}\text{ξ}\\ c^{-1}\text{ζ}& {\mu }_0\text{μ}\end{array}\right)$

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 $A=\frac{1}{\iota l}\ln \left(T\right)$, which is $A=\frac{1}{\iota l}U\left(\ln \left(P\right)\right)U^{-1}$ where $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 $k_{\text{a}}l$, $k_{\text{b}}l$, $k_{\text{c}}l$ or $k_{\text{d}}l$ is above $\pi$ or below $-\pi$, 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$\text{ε=}\left(\begin{array}{cc}\epsilon & 0\\ 0& \epsilon \end{array}\right)$,$\mu \text{=}\left(\begin{array}{cc}\mu & 0\\ 0& \mu \end{array}\right)$,$\text{ξ=}0$,$\zeta \text{=}0$surrounded by free space of impedance $z_0=\sqrt{{\mu }_0/{\epsilon }_0}$. 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:

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:

It’s clear that these eigenmodes correspond to forward and backward x- and y-polarized plane waves. Applying the formula in (20):

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

 

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.

Download Full Size | PDF

The diagonal terms of $\text{ξ}$ and $\zeta$ 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 y-polarized waves. Therefore, the model of this structure should contain only the off-diagonal terms of $\text{ξ}$ and $\zeta$ 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 $\text{ξ}$ and $\zeta$. 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.

 

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.

Download Full Size | PDF

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. $S_{\text{12}}^{\text{xx}}=n_2{S}_{\text{21}}^{\text{xx}}$ and $S_{\text{12}}^{\text{yy}}=n_2{S}_{\text{21}}^{\text{yy}}$, where $n_2$ is the refractive index of the backward substrate [18]), and this would result in having ${\zeta }_{\text{xy}}=-{\xi }_{\text{yx}}$ and ${\zeta }_{\text{yx}}=-{\xi }_{\text{xy}}$ as mathematically proved in [18], so we needed only to show ${\xi }_{\text{xy}}$ and ${\xi }_{\text{yx}}$ 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 $\text{ξ}$and $\zeta$.This structure is exactly the same from both sides, and this directional symmetry causes the off-diagonal terms of $\text{ξ}$and $\zeta$to be zero.

 

Fig. 4 Retrieval results of the two rod structure.

Download Full Size | PDF

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