Efficient Analysis of Metasurfaces in Terms of Spectral-Domain GSTC Integral Equations

TL;DR

This paper introduces a spectral-domain method for efficiently analyzing metasurfaces using GSTCs, reducing computational complexity by focusing on the surface rather than the entire domain, and demonstrating its effectiveness with practical examples.

Contribution

The paper develops a spectral-domain integral equation approach for metasurface analysis that is more efficient than traditional volumetric methods, focusing on surface equations.

Findings

Significant reduction in memory and computation time compared to finite-difference and finite-element methods.

Successful application to generalized-refractive and focusing metasurfaces.

Validation of the method's accuracy with illustrative examples.

Abstract

We present a spectral-domain (SD) technique for the efficient analysis of metasurfaces. The metasurface is modeled by generalized sheet transition conditions (GSTCs) as a zero-thickness sheet creating a discontinuity in the electromagnetic field. The SD expression of these GSTCs for a specified incident field leads to a system of four surface integral equations for the reflected and transmitted fields, which are solved using the method of moments in the spectral domain. Compared to the finite-difference and finite-element techniques that require meshing the entire computational domain, the proposed technique reduces the problem to the surface of the metasurface, hence eliminating one dimension and providing substantial benefits in terms of memory and speed. A monochromatic generalized-refractive metasurface and a polychromatic focusing metasurface are presented as illustrative examples.

Equations

$$\hat{\mathbf{z}}\times \Delta \mathbf{H}\left(\mathbf{\rho },t\right)$$

$$\hat{\mathbf{z}}\times \Delta \mathbf{E}\left(\mathbf{\rho },t\right)$$

$$\Delta \mathbf{E}\left(\mathbf{\rho },t\right)$$

$$\Delta \mathbf{H}\left(\mathbf{\rho },t\right)$$

$$\mathbf{P}\left(\mathbf{\rho },\omega \right)=$$

$$+\sqrt{{\mu }_0{ϵ}_0}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}_{\text{em}}\left(\mathbf{\rho },\omega \right)\cdot {\mathbf{H}}_{\text{av}}\left(\mathbf{\rho },\omega \right),$$

$$\mathbf{M}\left(\mathbf{\rho },\omega \right)=$$

$$+{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}_{\text{mm}}\left(\mathbf{\rho },\omega \right)\cdot {\mathbf{H}}_{\text{av}}\left(\mathbf{\rho },\omega \right),$$

$${\mathbf{E}}_{\text{av}}\left(\mathbf{\rho },\omega \right)$$

$${\mathbf{H}}_{\text{av}}\left(\mathbf{\rho },\omega \right)$$

$$\mathbf{\Psi }\left(\mathbf{\rho },t\right)=$$

$$\mathbf{\Psi }\left(\mathbf{\rho },\omega \right)=$$

$${\mathbf{E}}^{\text{i}}\left(\mathbf{r},t\right)$$

$${\mathbf{E}}^{\text{r}}\left(\mathbf{r},t\right)$$

$${\mathbf{E}}^{\text{t}}\left(\mathbf{r},t\right)$$

$$\int \cdots d\omega d{\mathbf{k}}_T={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }\cdots d\omega d{k}_{x}d{k}_y.$$

$$\nabla$$

$${\nabla }_T$$

$$k_z=\sqrt{k_0^2-k_x^2-k_y^2},$$

$${\mathbf{H}}^{\text{i}}\left(\mathbf{r},t\right)=\frac{1}{\omega {\mu }_0}\int$$

$${\tilde{\mathbf{E}}}^{\text{i}}\left({\mathbf{k}}_T,\omega \right)e^{j(\omega t-{\mathbf{k}}_T\cdot \mathbf{\rho }-k_{z}z)}d\omega d{\mathbf{k}}_T,$$

$${\mathbf{H}}^{\text{r}}\left(\mathbf{r},t\right)=\frac{1}{\omega {\mu }_0}\int$$

$${\tilde{\mathbf{E}}}^{\text{r}}\left({\mathbf{k}}_T,\omega \right)e^{(j\omega t-{\mathbf{k}}_T\cdot \mathbf{\rho }+k_{z}z)}d\omega d{\mathbf{k}}_T,$$

$${\mathbf{H}}^{\text{t}}\left(\mathbf{r},t\right)=\frac{1}{\omega {\mu }_0}\int$$

$${\tilde{\mathbf{E}}}^{\text{t}}\left({\mathbf{k}}_T,\omega \right)e^{j(\omega t-{\mathbf{k}}_T\cdot \mathbf{\rho }-k_{z}z)}d\omega d{\mathbf{k}}_T.$$

$$\mathbf{P}\left(\mathbf{\rho },\omega \right)$$

$$\mathbf{M}\left(\mathbf{\rho },\omega \right)$$

$${\overline{\stackrel{ˉ}{\mathbf{\chi }}}}_{\text{ee}}\left(\mathbf{\rho },t\right)$$

$${\overline{\stackrel{ˉ}{\mathbf{\chi }}}}_{\text{mm}}\left(\mathbf{\rho },t\right)$$

$${\tilde{\mathbf{P}}}_T\left({\mathbf{k}}_T,\omega \right)={\epsilon }_0\int {\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{ee}}^T\left({\mathbf{k}}_T^{\prime },\omega \right){\tilde{\mathbf{E}}}_{\text{av}}\left({\mathbf{k}}_T-{\mathbf{k}}_T^{\prime },\omega \right)d{\mathbf{k}}_T^{\prime },$$

$${\tilde{\mathbf{M}}}_T\left({\mathbf{k}}_T,\omega \right)=\int {\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{mm}}^T\left({\mathbf{k}}_T^{\prime },\omega \right){\tilde{\mathbf{H}}}_{\text{av}}\left({\mathbf{k}}_T-{\mathbf{k}}_T^{\prime },\omega \right)d{\mathbf{k}}_T^{\prime }.$$

$${\tilde{P}}_z\left({\mathbf{k}}_T,\omega \right)={\epsilon }_0\int {\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{ee}}^z\left({\mathbf{k}}_T^{\prime },\omega \right){\tilde{\mathbf{E}}}_{\text{av}}\left({\mathbf{k}}_T-{\mathbf{k}}_T^{\prime },\omega \right)d{\mathbf{k}}_T^{\prime },$$

$${\tilde{M}}_z\left({\mathbf{k}}_T,\omega \right)=\int {\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{mm}}^z\left({\mathbf{k}}_T^{\prime },\omega \right){\tilde{\mathbf{H}}}_{\text{av}}\left({\mathbf{k}}_T-{\mathbf{k}}_T^{\prime },\omega \right)d{\mathbf{k}}_T^{\prime }.$$

$${\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{ee}}^T\left({\mathbf{k}}_T,\omega \right)$$

$$\displaystyle\left[\begin{array}[]{ccc}\chi_{ee}^{xx}\left(\mathbf{k}_{T},\omega\right)&\chi_{ee}^{xy}\left(\mathbf{k}_{T},\omega\right)&\chi_{ee}^{xz}\left(\mathbf{k}_{T},\omega\right)\\ \chi_{ee}^{yx}\left(\mathbf{k}_{T},\omega\right)&\chi_{ee}^{yy}\left(\mathbf{k}_{T},\omega\right)&\chi_{ee}^{yz}\left(\mathbf{k}_{T},\omega\right)\end{array}\right],$$

$${\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{mm}}^T\left({\mathbf{k}}_T,\omega \right)$$

$$\displaystyle\left[\begin{array}[]{ccc}\chi_{\text{mm}}^{xx}\left(\mathbf{k}_{T},\omega\right)&\chi_{\text{mm}}^{xy}\left(\mathbf{k}_{T},\omega\right)&\chi_{\text{mm}}^{xz}\left(\mathbf{k}_{T},\omega\right)\\ \chi_{\text{mm}}^{yx}\left(\mathbf{k}_{T},\omega\right)&\chi_{\text{mm}}^{yy}\left(\mathbf{k}_{T},\omega\right)&\chi_{\text{mm}}^{yz}\left(\mathbf{k}_{T},\omega\right)\end{array}\right],$$

$${\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{ee}}^z\left({\mathbf{k}}_T,\omega \right)$$

$$\displaystyle\left[\begin{array}[]{ccc}\chi_{ee}^{zx}\left(\mathbf{k}_{T},\omega\right)&\chi_{ee}^{zy}\left(\mathbf{k}_{T},\omega\right)&\chi_{ee}^{zz}\left(\mathbf{k}_{T},\omega\right)\end{array}\right],$$

$${\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{mm}}^z\left({\mathbf{k}}_T,\omega \right)$$

$$\displaystyle\left[\begin{array}[]{ccc}\chi_{\text{mm}}^{zx}\left(\mathbf{k}_{T},\omega\right)&\chi_{\text{mm}}^{zy}\left(\mathbf{k}_{T},\omega\right)&\chi_{\text{mm}}^{zz}\left(\mathbf{k}_{T},\omega\right)\end{array}\right].$$

$$\mathbf{z}\times {\nabla }_T{P}_z\left(\mathbf{\rho },t\right)=\int \hat{\mathbf{z}}\times {\mathbf{k}}_T{\tilde{P}}_z\left({\mathbf{k}}_T,\omega \right)e^{j(\omega t-{\mathbf{k}}_T\cdot \mathbf{\rho })}d\omega d{\mathbf{k}}_T,$$

$$\mathbf{z}\times {\nabla }_T{M}_z\left(\mathbf{\rho },t\right)=\int \hat{\mathbf{z}}\times {\mathbf{k}}_T{\tilde{M}}_z\left({\mathbf{k}}_T,\omega \right)e^{j(\omega t-{\mathbf{k}}_T\cdot \mathbf{\rho })}d\omega d{\mathbf{k}}_T.$$

$$\begin{array}{rl}& \frac{1}{\omega {\mu }_0}\hat{\mathbf{z}}\times \left(-\left({\mathbf{k}}_T+k_z\hat{\mathbf{z}}\right)\times {\tilde{\mathbf{E}}}^{\text{i}}\left({\mathbf{k}}_T,\omega \right)-\left({\mathbf{k}}_T-k_z\hat{\mathbf{z}}\right)\times {\tilde{\mathbf{E}}}^{\text{r}}\left({\mathbf{k}}_T,\omega \right)+\left({\mathbf{k}}_T+k_z\hat{\mathbf{z}}\right)\times {\tilde{\mathbf{E}}}^{\text{t}}\left({\mathbf{k}}_T,\omega \right)\right)=\\ & \frac{1}{2}i\omega {\epsilon }_0\int {\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{ee}}^T\left({\mathbf{k}}_T-{\mathbf{k}}_T^{\prime },\omega \right)\left[{\tilde{\mathbf{E}}}^{\text{i}}\left({\mathbf{k}}_T^{\prime },\omega \right)+{\tilde{\mathbf{E}}}^{\text{r}}\left({\mathbf{k}}_T^{\prime },\omega \right)+{\tilde{\mathbf{E}}}^{\text{t}}\left({\mathbf{k}}_T^{\prime },\omega \right)\right]d{\mathbf{k}}_T^{\prime }\\ & -\frac{1}{2\omega {\mu }_0}\int \hat{\mathbf{z}}\times \left({\mathbf{k}}_T^{\prime }-{\mathbf{k}}_T\right){\stackrel{~}{\overline{\stackrel{ˉ}{\mathbf{\chi }}}}}_{\text{mm}}^z\left({\mathbf{k}}_T^{\prime }-{\mathbf{k}}_T,\omega \right)\left[\left({\mathbf{k}}_T^{\prime }+k_z^{\prime }\hat{\mathbf{z}}\right)\times \left[{\tilde{\mathbf{E}}}^{\text{i}}\left({\mathbf{k}}_T^{\prime },\omega \right)+{\tilde{\mathbf{E}}}^{\text{t}}\left({\mathbf{k}}_T^{\prime },\omega \right)\right]+\left({\mathbf{k}}_T^{\prime }-k_z^{\prime }\hat{\mathbf{z}}\right)\times {\tilde{\mathbf{E}}}^{\text{r}}\left({\mathbf{k}}_T^{\prime },\omega \right)\right]d{\mathbf{k}}_T^{\prime },\end{array}$$

Full text

I Introduction

Metasurfaces are 2D arrangements of scattering particles that transform electromagnetic waves in specified ways. In contrast to 3D metamaterials, they are easier to fabricate, and compared to frequency selective surfaces (FSSs), they provide richer transformations. Metasurfaces have already found a great diversity of applications, including polarization transformers and rotators [1, 2], perfect absorbers [3], nonreciprocal screens [2, 4, 5, 6], beam transformers [7], electromagnetic interferometers and processors [8].

A metasurface can be effectively modeled as a zero-thickness sheet that creates a discontinuity in the electromagnetic field. Such a discontinuity can be most generally described by generalized sheet transition conditions (GSTCs) [9, 10, 11, 12], where the discontinuity in the electric and magnetic fields are related to surface electric and magnetic current and surface electric and magnetic polarization densities. GSTCs are not implemented in currently available commercial software and their future implementation will require efficient numerical approaches. Exact zero-thickness modeling by GSTCs has been recently incorporated in academic numerical solvers based on finite differences [13], including finite-difference time-domain (FDTD) [14], finite-difference frequency-domain (FDFD) [15], and finite-element methods [16, 17]. The FDTD-GSTC and FDFD-GSTC techniques [13] properly model the zero thickness of the metasurface through the introduction of virtual nodes at the location of the metasurface on, while the FEM-GSTC does this naturally through flexible FEM mesh points at both sides of the metasurface.

Although they model the zero thickness electromagnetic field discontinuity quite well, the aforementioned GSTC finite-difference and finite-element techniques require to mesh the entire computational domain. In large problems, and particularly three-dimensional ones, volumetric meshes represent a high computational burden in terms of memory and speed. In contrast, surface integral equations reduce the problem complexity by one dimension, and are hence able to solve large problems that may be intractable with other techniques. As GSTCs describe the electromagnetic discontinuity on a surface, it is natural to cast them to the form of surface integral equations [18, 19] and thus essentially reduce the problem to the surface of the metasurface. Moreover, for a flat metasurface, GSTCs may be conveniently expressed in the spectral domain (SD) [20] and analyzed with the well developed Fourier techniques [20]. This paper combines these two points to provide an exact analysis of metasurfaces in the spectral domain (SD) [20] based on SD integral equations. GSTCs are modeled as exact zero-thickness transition conditions and transformed into a system of SD surface integral equations which are then solved using the method of moments (MoM) [18]. The integral equations are derived and solved in the temporal frequency domain. However, the method can be generalized to time domain SD integral equation as well if the problem is linear.

The organization of the paper is as follows. GSTCs are described in Sec. II. SD integral equations are derived in Sec. III. MoM solution of the SD integral equations are presented in Sec. IV. Two illustrative examples are provided in Sec. V. Finally, Sec. VI concludes the paper.

II Transition Conditions

Figure 1 represents a general metasurface transformation, where a given incident wave is transformed into specified reflected and transmitted waves. Assuming that the metasurface is placed in the $xy$-plane at $z=0$, the electromagnetic field discontinuity it induces may be rigorously modeled by the GSTCs [9, 11]

[TABLE]

where $\mathbf{P}$ and $\mathbf{M}$ represent the generally space-time dependent electric and magnetic surface polarization densities, and where

[TABLE]

with superscripts + and - representing the total fields on the right and left of the metasurface, respectively, subscripts T and z representing the components tangential and normal to metasurface, respectively, and $\boldsymbol{\rho}$ representing an arbitrary point on the metasurface, i.e. $\boldsymbol{\rho}=x\hat{\mathbf{x}}+y\hat{\mathbf{y}}$.

Assuming that the metasurface is dispersive but linear and local, the polarization densities are related to the average (${}_{\text{av}}$) electric and magnetic fields on the metasurface through the general bianisotropic dispersion relations

[TABLE]

with

[TABLE]

where $\omega$ represents the temporal frequency. The functions of $\left(\boldsymbol{\rho},\omega\right)$ in these relations are related to their temporal Fourier transforms $\left(\boldsymbol{\rho},t\right)$ functions by the Fourier transform pair

[TABLE]

where $\boldsymbol{\Psi}=\mathbf{P}$, $\mathbf{M}$, $\mathbf{E}$ and $\mathbf{H}$. The next section will derive the GSTCs in the spatiotemporal $(\mathbf{k},\omega)$ spectral domain and develop the corresponding SD integral equations for the reflected and transmitted fields.

III Spectral-Domain Analysis

It is assumed that incident field and metasurface susceptibilities are known. This section derives the SD integral equations for the corresponding reflected and transmitted fields. The known incident field and the unknown reflected and transmitted fields are decomposed into their spatiotemporal Fourier components $(\mathbf{k}_{T},\omega)$ as follows:

[TABLE]

where $\mathbf{k}_{T}=k_{x}\hat{\mathbf{x}}+k_{y}\hat{\mathbf{y}}$ and where the integral sign specifically represent the triple integral

[TABLE]

With these definitions, the differential operators $\nabla$ and $\nabla_{T}$ reduce to their algebraic SD equivalents

[TABLE]

where $k_{z}$ is found as

[TABLE]

with $k_{0}=\omega_{0}\sqrt{\epsilon_{0}\mu_{0}}$, upon inserting (6) into the free-space Helmholtz equation, $(\nabla^{2}+k_{0}^{2})\mathbf{E}\left(\mathbf{r},\omega\right)=0$, and decomposing the result $k^{2}+k_{0}^{2}=0$ into transverse and normal components. The sign of the square root must be chosen so as to satisfy the radiation condition at $z=\pm\infty$, i.e. $\Re(k_{z})\geq 0$, $\Im(k_{z})\leq 0$.

The magnetic fields corresponding to (6) are obtained through the Maxwell-Faraday equation, $\mathbf{H}\left(\mathbf{r},\omega\right)=\frac{-j\omega}{\mu_{0}}\nabla\times\mathbf{E}\left(\mathbf{r},\omega\right)$, as

[TABLE]

We next assume, for simplicity and without loss of generality, monoanisotropic (rather than bianisotropic) metasurfaces, i.e.

[TABLE]

where the electric and magnetic susceptibility tensors are also expanded in the spectral domain:

[TABLE]

In the spatial-frequency Fourier domain the product expressions in (11) become convolutions, and the corresponding transverse and normal polarization densities read thus

[TABLE]

and

[TABLE]

where

[TABLE]

Finally, the terms involving $\hat{\mathbf{z}}\times\nabla_{T}$ in (1) become in SD form

[TABLE]

At this point, substituting (6), (10), (13) and (16) into (1), and using the orthogonality property of Fourier harmonics, yields the following SD integral equations:

Equations (17) and (18) represent four integral equations in six unknowns $\tilde{E}^{\text{r}}_{x}$, $\tilde{E}^{\text{r}}_{y}$, $\tilde{E}^{\text{r}}_{z}$, $\tilde{E}^{\text{t}}_{x}$, $\tilde{E}^{\text{t}}_{y}$, $\tilde{E}^{\text{t}}_{z}$. However, as $\mathbf{E}^{\text{r}}$ and $\mathbf{E}^{\text{t}}$ satisfy the free-space Maxwell equations, the normal components $\tilde{E}^{\text{r}}_{z}$ and $\tilde{E}^{\text{t}}_{z}$ are related to the transverse components through the divergence relations $\nabla\cdot\mathbf{E}^{\text{r}}(\mathbf{r},t)=0$, $\nabla\cdot\mathbf{E}^{\text{t}}(\mathbf{r},t)=0$ and are thus eliminated using

[TABLE]

Substituting (19) into (17) and (18) finally results in a system of four coupled integral equations in the four unknowns $\tilde{E}^{\text{r}}_{x}(\mathbf{k}_{T},\omega)$, $\tilde{E}^{\text{r}}_{y}(\mathbf{k}_{T},\omega)$, $\tilde{E}^{\text{t}}_{x}(\mathbf{k}_{T},\omega)$, $\tilde{E}^{\text{t}}_{y}(\mathbf{k}_{T},\omega)$, that may be solved by the method of moments.

IV Integral-Equation System Resolution by the Method of Moments

The SD integral-equation system (17) and (18) with (19) are expressed in the temporal frequency ($\omega$) domain. However, if the system is linear, they naturally also apply to non time harmonic problems, since any (non time harmonic) field may be analyzed as the superposition of its time-harmonic Fourier components. In particular, a pulse incident field may be written as

[TABLE]

For this excitation, we first transform $\mathbf{E}^{\text{i}}\left(\boldsymbol{\rho},\omega\right)$, $\bar{\bar{\boldsymbol{\chi}}}_{\text{ee}}\left(\boldsymbol{\rho},\omega\right)$ and $\bar{\bar{\boldsymbol{\chi}}}_{\text{mm}}\left(\boldsymbol{\rho},\omega\right)$ into their spatial Fourier harmonics, namely $\tilde{\mathbf{E}}^{\text{i}}\left(\mathbf{k}_{T},\omega\right)$, $\tilde{\bar{\bar{\boldsymbol{\chi}}}}_{\text{ee}}\left(\mathbf{k}_{T},\omega\right)$ and $\tilde{\bar{\bar{\boldsymbol{\chi}}}}_{\text{mm}}\left(\mathbf{k}_{T},\omega\right)$, which are used as (known) inputs to the integral equation system. The system is then solved so as to determine the unknowns $\tilde{E}^{\text{r}}_{x}(\mathbf{k}_{T},\omega)$, $\tilde{E}^{\text{r}}_{y}(\mathbf{k}_{T},\omega)$, $\tilde{E}^{\text{t}}_{x}(\mathbf{k}_{T},\omega)$ and $\tilde{E}^{\text{t}}_{y}(\mathbf{k}_{T},\omega)$. The process is repeated for all the frequencies $\omega$ of the input wave111No new frequencies are created since the system is assumed to be both linear and time-invariant. spectrum. Finally, the space-time dependent reflected and transmitted fields are found using (6). In the case of a time-harmonic incident wave, the system is naturally solved only once, for the incident frequency.

For each frequency sample $\omega$, the SD integral-equation system (17) and (18) with (19) is solved using a standard method of moment (MoM) technique in the spectral $\mathbf{k}_{T}$ domain [18, 21, 22]. First, the unknowns are expanded over a set of basis functions, $B_{n}(\mathbf{k}_{T})$, as

[TABLE]

where $a_{n}^{\text{r}x}$, $a_{n}^{\text{r}y}$, $a_{n}^{\text{t}x}$, $a_{n}^{\text{t}y}$ represent unknown coefficients. Equations (17) and (18) are then multiplied by a set of SD testing functions $W_{n}(\mathbf{k}_{T})$ and integrated over the entire $\mathbf{k}_{T}$ domain. This process reduces (17) and (18) with (19) to a system of linear equations for the unknown coefficients $a_{n}^{\text{r}x}$, $a_{n}^{\text{r}y}$, $a_{n}^{\text{t}x}$, $a_{n}^{\text{t}y}$. In the forthcoming illustrative examples, we have used point matching, i.e. $W_{n}(\mathbf{k}_{T})=\delta(\mathbf{k}_{T}-\mathbf{k}_{nT})$, and Galerkin testing, i.e. $B_{n}(\mathbf{k}_{T})=W_{n}(\mathbf{k}_{T})$.

V Illustrative Examples

This section presents two illustrative examples of application of the proposed SD integral-equation resolution technique. Section V-A considers a monochromatic generalized-refractive metasurface, while Sec. V-B considers a polychromatic focusing metasurface, with different focal points at different frequencies. In both cases, we first obtain the metasurface susceptibility functions for the specified transformation using the GSTC synthesis technique presented in [10], and next analyze the resulting metasurface using the SD integral-equation technique and verify that the synthesized operation is properly achieved.

V-A Monochromatic Generalized-Refractivie Metasurface

The monochromatic generalized-refractive metasurface represented in Fig. 4. It transforms a given incident plane wave with angle $\theta^{\text{i}}$ into the specified reflected and transmitted waves, with transmission and reflection coefficients $T$ and $R$ and angles $\theta^{\text{r}}$ and $\theta^{\text{t}}$, respectively. Assuming that the metasurface is monoisotropic and devoid of normal susceptibilities, its susceptibility tensors reduce to

[TABLE]

corresponding to the tangential and normal polarization densities

[TABLE]

We require that the incident wave impinge on the metasurface with the angle $\theta^{i}=60^{\circ}$, negatively refract with the angle $\theta^{t}=60^{\circ}$ and transmission coefficient $T=0.8$, and reflect with the angle $\theta^{r}=25\deg$ and reflection coefficient $R=-0.6$. The metasurface satisfying to these specifications is readily synthesized from the GSTCs relations

[TABLE]

where the difference and average fields are composed of the specified fields according to (2) and (4). Since the specified fields have no $x$ dependence, the susceptibilities are only functions of $y$. They are plotted in Fig. 3. Note that as this transformation conserves the total power $|T|^{2}+|R|^{2}=1$, losses in the susceptibilities should be compensated by gains. Therefore in some region the imaginary part of susceptibilities become positive, corresponding to gain.

From this point, we shall carry out the analysis and see if we thereby retrieve the synthesized operation. Given the incident field and synthesized metasurface susceptibilities, we compute the reflected and transmitted fields using the SD integral-equation and MoM technique. The results are plotted in Fig. 4. They exactly correspond to the specified fields, which validates the analysis technique.

V-B Polychromatic Focusing Metasurface

The polychromatic focusing metasurface is designed so as to exhibit different focal points at different frequencies, as shown in Fig. 5. A plane wave multi-wavelength pulse normally impinges on the metasurface from the left, and the metasurfaces focuses the different wavelengths or frequencies onto corresponding different points on the right. Such a metasurface essentially operates as a real-time spectrum analyzer as it spatially decomposes the frequency contents of the input pulse in real time. To synthesize the metasurface, we make separate designs for each of the frequencies $\omega_{i}$ and calculate the corresponding susceptibilities $\chi_{\text{ee}}\left(y,\omega_{i}\right)$, $\chi_{\text{mm}}\left(y,\omega_{i}\right)$ from the specified electromagnetic fields at the left and right sides of the metasurface.

Consider the specified field on the left side of the metasurface to be a $y$-polarized normally incident plane wave with frequency $\omega_{i}$ designed to create a $y$-polarized focus at the corresponding specified focal point. The required field on the right side of the metasurface can be calculated by placing a $y$-polarized dipole at the frequency $\omega_{i}$ on the focal point, calculating the field radiated by this dipole on the metasurface at $z=0_{+}$ and applying phase conjugation to this field to achieve reverse propagation. The susceptibilities $\chi_{\text{ee}}\left(y,\omega_{i}\right)$, $\chi_{\text{mm}}\left(y,\omega_{i}\right)$ are then computed as in Sec. V-A. Again, as in (23), we assume purely transverse and monoisotropic susceptibilities.

The so obtained susceptibilities are plotted in Fig. 6 for a four-wavelength metasurface with equidistant frequencies across the range $\omega_{1}\cdots 2\omega_{1}$ and corresponding four focal points. The horizontal axis is normalized to the wavelength at frequency $\omega_{1}$. The focal points are separated by $1.33\lambda_{1}$ from $-2\lambda_{1}$ to $2\lambda_{1}$. Note that the susceptibilities exhibit maxima at the points corresponding to their respective focal height. The electric and magnetic susceptibilities are almost identical, corresponding to a matched, and hence reflection-less, metasurface [23].

Next we perform the SD integral-equation and MoM analysis of the metasurface just synthesized. The sum of the resulting reflected and transmitted fields at frequencies $\omega_{1}$ to $\omega_{4}$ are plotted in Fig. 7. The metasurface is well matched (no reflection) and generates focal spots at the four specified points. Moreover, these spots correspond to the specified frequencies as may deduced from the corresponding wavelengths. The focal spots do not seem very accurate in Fig. 7, but this is due to the superposition of the four waves of frequencies $\omega_{1}$ to $\omega_{4}$ and this does not represent a practical issue since the different frequencies are orthogonal and hence independent from each other. Figure 8, that shows the response of the metasurface separately excited by the four waves, confirms this fact.

VI Conclusion

We have presented a spectral domain surface integral equation technique for the analysis of metasurfaces. A system of four surface integral equations have been derived for the reflected and transmitted electric fields, and solved using the method of moments. Compared to the finite difference and finite element techniques that require volumetric meshes, the proposed technique reduces the problem to the surface of the metasurface, eliminating one dimension and therefore providing benefits in terms of memory and speed. A monochromatic generalized refraction metasurface and a polychromatic focusing metasurface have been presented as illustrative examples.