Ray-Optical Evaluation of Scattering from Electrically Large Metasurfaces Characterized by Locally Periodic Surface Susceptibilities

Electromagnetic (EM) metasurfaces have recently emerged as spatial field processors, which can be engineered to transform incident wavefronts into desired wavefronts in both reflection and transmission. They are thin surfaces (with thicknesses much smaller than the wavelength, $\lambda$) consisting of an array of sub-wavelength scatterers (or particles) which interact with the incident waves, and through their rich electric and magnetic polarization responses, can manipulate their amplitude, phase and polarizations. Consequently, they have found a myriad of applications across the entire electromagnetic spectrum ranging from optics to microwaves [2, 3, 4, 5, 6].

While they are created using engineered sub-wavelength unit-cell structures, practical metasurfaces are electrically large – of the order of several tens of wavelengths or more. The incident field interacts with the non-uniformly arrayed resonating particles to generate a spatial distribution of localized microscopic responses which together with the edge effects due to the finite size of the surface aperture, produce the overall macroscopic response. In addition, if the metasurface is part of a larger system of scattering objects, the field interactions become even more complex and the scattering problem must be seen as a multi-scale field problem – i.e., from the sub-wavelength unit-cell level to the system level.

The implications of the multi-scale nature of the metasurface scattering problem may be better understood by considering the application of metasurfaces to optimizing wireless communication performance at microwave and mmWave frequencies. For instance, radio-frequency (RF) metasurfaces are envisioned to act as smart scatterers manipulating EM wave propagation in a given radio environment, avoiding obstacles or redirecting wireless signals to desired receiver locations using the best possible channel path [7, 8, 9]. Another application to consider is that of EM illusions and camouflaging [10, 11], where electrically large metasurfaces can act as holograms to project a virtual image of an object, or to enclose a real object, camouflaging it against a complex scattering environment. In all cases, such environments are naturally electrically large, with a variety of static and moving objects, of sizes which are orders of magnitude larger than those of the metasurface unit cells.

The problem of determining the scattered fields from an electrically large surface, possibly in the presence of other scattering objects, is a challenging task. Among the variety of numerical methods that have recently been proposed to model field scattering from metasurfaces are FDTD and FDFD [12, 13, 14, 15, 16], FEM [17], Spectral Domain method [18] and Integral Equation Solvers (IE-GSTCs) based on Boundary Element Methods [19]. Since volumetric numerical methods (FDFD, FDTD and FEM) based on finite domains are too inefficient to model electrically large metasurfaces, surface based methods such as IE-GSTCs become a preferred choice. In the surface based approach, metasurface structures are conveniently represented by equivalent zero-thickness models to reduce computational complexity while preserving their physical scattering characteristics. The EM properties of the zero-thickness sheet can then be described in terms of surface susceptibility tensors $\bar{\bm{\chi}}$ which, in conjunction with the Generalized Sheet Transition Conditions (GSTCs), provide a rigorous field model from which to determine the scattered surface fields anywhere in space [20, 21, 22].

While surface based methods have clear computational complexity advantages, for very large metasurface problems (or alternatively, very high frequencies) they also become unworkable due to excessive resource requirements. If concessions can be made on the desired level of accuracy, ray methods based on geometrical optics (GO) and a suitable uniform asymptotic description of diffraction effects, such as the uniform theory of diffraction (UTD) [23], emerge as a promising alternative. UTD based (or similar) high-frequency methods are computationally efficient and present an intuitive picture of the fields scattered from the surface. Consequently, they have been used in a wide variety of applications including radio propagation modeling [24] and analysis of EM structures such as impenetrable gratings and diffusers [25], conventional curved surfaces [26], and diffraction gratings with slowly varying parameters [27].

Motivated by such applications, a uniform asymptotic description of physical optics (PO) surface scattering has recently been applied to metasurfaces [1]. In [1], a generic metasurface is modeled as a zero-thickness sheet with prescribed reflection and transmission coefficients. Assuming that the metasurface is locally periodic, these coefficients and the corresponding scattered surface fields are expanded into their spatial Fourier components, which are then incorporated in the EM field integral equations. Evaluation of these field integrals in the high-frequency limit ($k\rightarrow\infty$, $k=2\pi/\lambda$ being the wavenumber in free space) leads to an efficient computation of the scattered fields around the metasurface using contributions from a discrete set of critical points. While the ray-based numerical framework is computationally efficient, the prescribed reflection and transmission coefficients are in reality field-dependent, and should therefore be respecified or otherwise updated every time the incident field direction is changed.

It is more desirable to define the scattering properties of the metasurface in terms of its intrinsic parameters, such as surface susceptibilities $\bar{\bm{\chi}}$, from which the fields can be computed on demand, for an arbitrary incidence angle. This work, therefore, represents a continuation of the uniform ray description of physical optics scattering by metasurfaces presented in [1], in which the GSTCs are directly incorporated into the ray formulation, and the metasurface is now described in terms of angle-independent surface susceptibilities $\bar{\bm{\chi}}$. This approach combines the elegance of the numerical framework, the local periodicity of the metasurface and rigorous GSTC boundary conditions to model the wave transformation capabilities of metasurfaces based on surface susceptibilities.

The paper is organized as follows: Sec. II describes the problem statement and the GSTC model of a generic metasurface. Sec. III then modifies the GSTCs for the case of GO fields interacting with a locally periodic metasurface in the high-frequency limit. This entails (locally) expressing the prescribed surface susceptibility profile in terms of a Fourier series expansion, and the scattered field as plane-wave spectra of spatial modes. Equations are provided for numerically solving the scattered surface fields for a given GO incident field, and propagating them to arbitrary detector locations away from the surface. A simple procedure to convert an arbitrarily prescribed susceptibility profile into the desired Fourier series form is described in Sec. IV. This procedure uses a spatially windowed Fourier transform, and is similar to the conventional use of Short Time Fourier Transforms (STFT). Several numerical examples are then presented in Sec. V to demonstrate the proposed ray-optical GSTC (RO-GSTC) formulation. Conclusions are finally provided in Sec. VI.

Consider a planar, non-uniform metasurface in the plane spanned by the unit basis vectors $\hat{\bf x}$ and $\hat{\bf y}$, whose normal vector, $\hat{\bf n}$, is identical to $-\hat{\bf z}$, as shown in Fig. 1. A known incident field excites the metasurface, and generates the scattered reflection and transmission fields in free-space Regions I and II, respectively. The properties of the metasurface are assumed to vary along $\hat{\bf x}$ only for simplicity, and the incident field is assumed to be cylindrical and in the $x$–$z$ plane.

The main objective of this paper is to compute the unknown amplitudes of the reflected and transmitted fields at all surface points $\mathbf{q}^{\prime}$, and then to evaluate the total scattered fields away from the metasurface as a superposition of field contributions originating from the surface points. The scattering response on the surface can be described in terms of localized transmission and reflection dyadics $\bar{\bm{\Gamma}}$ and $\bar{\bf T}$ each point on the surface [1], well-known to be dependent on the propagation direction of the incident field. Alternatively, its scattering properties can be described by a well-chosen set of constitutive parameters which are, by definition, field-independent. The constitutive parameters of choice for metasurfaces are the electric and magnetic surface susceptibility tensors $\bar{\bm{\chi}}_{\rm ee}$, $\bar{\bm{\chi}}_{\rm mm}$, $\bar{\bm{\chi}}_{\rm em}$ and $\bar{\bm{\chi}}_{\rm me}$, defined by up to 36 complex scalars that capture the microscopic response of the surface to the average fields around it.

For later convenience, we denote the perpendicular and parallel components of an arbitrary vector ${\bf A}$ with respect to the plane with normal $\hat{\bf n}$ by

where $\bar{\bf I}$ is the unit dyadic. The unit vectors $\hat{\bf u}_{1}$, $\hat{\bf u}_{2}$ and $\hat{\bf u}_{3}$ are defined to be identical to $\hat{\bf x}$, $\hat{\bf y}$ and $\hat{\bf z}$, respectively, and the components of a vector ${\bf A}$ along $\hat{\bf u}_{p}$, $p=\{1,2,3\}$, are denoted by

Similarly, the components of a dyadic $\bar{\bf A}$ are denoted by

Finally, the superscript $(m)$ denotes the m${}^{\text{th}}$ Fourier coefficient of any scalar, vector or dyadic quantity that is a periodic function of $\psi$ with period $\lambda$, i.e.,

As shown in Sec. III-C, the surface fields (8) needed to compute the scattered fields away from the metasurface can be determined under the assumption that the surface susceptibilities are periodic functions (4) of the parameter $\psi$. In (4), $\bar{\bm{\chi}}^{(m)}_{ab}(x^{\prime})$ denotes the amplitude of the $m^{th}$ mode of $\bar{\bm{\chi}}_{ab}$ near $x^{\prime}$, which has a localized spatial frequency $km\dot{\psi}(x^{\prime})$. Our objective now is to estimate these location-dependent spatial frequencies across an arbitrary given surface susceptibility profile, and determine the corresponding modal amplitudes. To this end, we investigate the Short Time Fourier Transform (STFT), conventionally applied to time-varying non-stationary signals [36], as a means to estimate the localized spatial frequency content of $\bar{\bm{\chi}}_{ab}$. When applied to the problem considered herein, we analogously refer to this procedure as the Short Space Fourier Transform (SSFT).

To describe the procedure, assume for the moment that the metasurface is characterized by a scalar susceptibility that is locally periodic, and can therefore be expressed as

where $\psi(x)$ is well approximated by (5) if $x$ is sufficiently close to $x^{\prime}$. The SSFT of $\chi$, denoted by $\mathcal{S}\{\chi(x)\}=\mathcal{S}(\kappa,x^{\prime})$, involves computing the Fourier transform of $\chi$ multiplied by a space-limited window function $w(x-x^{\prime})$ centered at $x=x^{\prime}$, which slides along the $x$-axis; $\kappa$, in rad/m, denotes the spatial frequency. At each window location,

represents the localized spectral content of the profile. Substituting (26) and (5), it is not difficult to show that $\mathcal{S}\{\chi\}$ is well approximated by

provided that the window function is narrow enough for the quadratic term in (5) to be negligibly small for all $x$ where $w(x-x^{\prime})\neq 0$. In (28),

is the window function in the spatial Fourier domain.

The SSFT distribution of a locally periodic susceptibility profile typically looks like that of Fig. 3, which shows discrete bands of high amplitude that are uniformly spaced along $\kappa$, with a spacing that varies with the window location $x^{\prime}$. Recognizing that these bands correspond with the summation terms of (28), their spacing, $\Delta\kappa$, equals $k\dot{\psi}$ and the locations of their maxima with respect to $\kappa$, denoted by $\kappa_{m}$, can be used to estimate the local periodicity of the metasurface as follows:

Having estimated $\dot{\psi}$ for a sufficiently dense grid of surface points, the phase function $\psi(x)$ can be numerically constructed via

where the constant of integration can be chosen arbitrarily, and is set to zero for simplicity here.

Generalizing this procedure to a set of dyadic surface susceptibility profiles $\bar{\bm{\chi}}_{ab}(x)$, with $a\in\{{\rm e},{\rm m}\}$, the gradient (slope) function $\dot{\psi}(x)$ is first determined through (29) using the localized spatial frequencies estimated from a dyadic version of the SSFT integral (27). The phase function $\psi(x)$ is next computed using (30) using with a polynomial curve fitted to the extracted $\dot{\psi}$ data. Care must be taken that the assumption of linearly varying phase within the spatial limits of the window function $w(\cdot)$ is satisfied. Finally, each unknown $\bar{\bm{\chi}}_{ab}^{(m)}(x)$ may be obtained by direct numerical integration of $\bar{\bm{\chi}}_{ab}(x)$ (see (1)), using $\bar{\bm{\chi}}_{ab}(\psi;x^{\prime})\simeq\bar{\bm{\chi}}_{ab}(x)$ and the estimated $\psi(x)$ from (30),

where $\psi_{1,2}$ lie on either side of $\psi(x^{\prime})$ and span one full period such that $\psi_{2}-\psi_{1}=\lambda$. We can then use the relationship $d\psi=\dot{\psi}dx$ to change the integration to be over $x$ and we obtain,

with the limits given by $x_{1,2}=\psi^{-1}(\psi_{1,2}$).

It is conceptually straightforward to extend the procedure to the wider class of 2D (locally periodic) metasurfaces, whose susceptibility profiles vary along the direction of a gradient vector that is not necessarily parallel to $\hat{\bf x}$. To do this, a 2D spatial Fourier transform needs to be applied to $\bar{\bm{\chi}}_{ab}(x,y)$ multiplied by a 2D window function. The locations $(\kappa_{1,m},\kappa_{2,m})$ of the spectral maxima at a given surface location correspond to a 2D spatial frequency vector ${\bm{\kappa}}_{m}=\kappa_{1,m}\hat{\mathbf{u}}_{1}+\kappa_{2,m}\hat{\mathbf{u}}_{2}$, from which one can construct the phase function $\psi(x,y)$ and its gradient ${\bf g}_{\psi}(x,y)$.

While the above procedure is simple, it also has limitations. The SSFT integral (27) suffers from the well-known uncertainty principle, according to which its resolution (the minimum distance at which nearby space-frequency components can be resolved) is inversely related along $\kappa$ and $x$. Resolution in one domain can be traded off for increased resolution in the other by choosing a different window function $w(\cdot)$, but in some cases it may not be possible to find a window function that avoids any significant space-frequency overlap. However, as we shall see, the use of a fitted function to represent $\dot{\psi}$ alleviates these issues if discretion is taken in the retention of points in areas where the SSFT cannot discriminate the modal frequencies clearly.

The metasurface is assumed to be perfectly matched (${\bf E}_{r}={\bf 0}$) and transmit an electric and magnetic surface field ${\bf E}_{t}(x)=0.2j\hat{\bf y}$ and ${\bf H}_{t}(x)=-0.2j/\eta\hat{\bf x}$, where $\eta$ is the free space impedance. The metasurface therefore introduces a position dependent loss that accommodates the variation in $|\mathbf{E}_{i}|$ along the surface to produce a uniform transmitted field intensity. If we constrain the surface synthesis to an isotropic surface, described by the tangential components of $\bar{\bm{\chi}}_{\rm mm}$ only, the corresponding electric and magnetic surface susceptibilities are then obtained using:

where $\Delta\{\cdot\}$ is the field difference across the metasurface, and $\{\cdot\}_{\text{avg}}$ is the average fields [28]. For this example, the computed $\chi_{\text{ee}}(x)$ ($\chi_{\text{mm}}$ not shown for conciseness) is shown in Fig. 4(a).

In order to express the susceptibilities in the form of (4), we first compute the SSFT using MATLAB’s STFT function, the results of which are shown in Fig. 4(b). The window function $w(p)$ used was the discrete Kaiser window [38][39]. The surface was discretized with a spacing of $\lambda/400$ and the spatial extent of the window is approximately 0.20 m. As can be seen, the SSFT decomposes the surface susceptibility into distinct spatial modes for multiple positions along the surface. Using the $m=1$ mode, the spatial frequencies are then extracted, and used in (29) to compute $\dot{\psi}$, which is shown in Fig. 4(c). Notice that for this particular example, it is difficult to extract $\kappa_{m}$ near $x=0$ due to limited spatial and frequency resolution.