Scattering Field Solutions of Metasurfaces based on the Boundary Element Method (BEM) for Interconnected Regions

TL;DR

This paper introduces a Boundary Element Method (BEM) approach to compute electromagnetic scattering from interconnected regions with metasurfaces, validated through Gaussian beam and diffraction grating examples.

Contribution

It develops a BEM-based framework incorporating GSTCs for analyzing EM scattering from complex metasurfaces in interconnected regions.

Findings

Validated with semi-analytical method for Gaussian beam scattering.

Confirmed results with FDFD for diffraction grating.

Demonstrated effectiveness of BEM in metasurface scattering problems.

Abstract

A methodology for determining the scattered Electromagnetic (EM) fields present for interconnected regions with common metasurface boundaries is presented. The method uses a Boundary Element Method (BEM) formulation of the frequency domain version of Maxwell's equations - which expresses the fields present in a region due to surface currents on the boundaries. Metasurface boundaries are represented in terms of surface susceptibilities which when integrated with the Generalized Sheet Transition Conditions (GSTCs), gave rise to an equivalent configuration in terms of electric and magnetic currents. Such a representation is then naturally incorporated into the BEM methodology. Two examples are presented for EM scattering of a Gaussian beam to illustrate the proposed method. In the first example, metasurface is excited with a diverging Gaussian beam, and the scattered fields are validated using a semi-analytical method. Second example concerned with a non-uniform metasurface modeling a diffraction grating, whose results were confirmed with conventional Finite Difference Frequency Domain (FDFD) method.

Equations

$$-\hat{\mathbf{n}}\times {\mathbf{H}}_{s,1}={\mathbf{J}}_s,$$

$$\hat{\mathbf{n}}\times {\mathbf{E}}_{s,1}={\mathbf{K}}_s.$$

$$-\hat{\mathbf{n}}\times \left({\mathbf{E}}_{s,1}-{\mathbf{E}}_{s,2}\right)={\mathbf{K}}_s$$

$$\hat{\mathbf{n}}\times \left({\mathbf{H}}_{s,2}-{\mathbf{H}}_{s,1}\right)={\mathbf{J}}_s,$$

$$\hat{\mathbf{n}}\times \Delta \mathbf{H}=j\omega {\mathbf{P}}_{||}-\hat{\mathbf{n}}\times {\nabla }_{||}{\mathbf{M}}_n$$

$$\hat{\mathbf{n}}\times \Delta \mathbf{E}=-j\omega {\mu }_0{\mathbf{M}}_{||}-\hat{\mathbf{n}}\times {\nabla }_{||}\left(\frac{{\mathbf{P}}_n}{{ϵ}_0}\right)$$

$$\mathbf{P}(\omega )=ϵ{\overline{\overline{\chi }}}_{\text{ee}}{\mathbf{E}}_{\text{avg}}(\omega )+\sqrt{\mu ϵ}{\overline{\overline{\chi }}}_{\text{em}}{\mathbf{H}}_{\text{avg}}(\omega )$$

$$\mathbf{M}(\omega )=\sqrt{\frac{ϵ}{\mu }}{\overline{\overline{\chi }}}_{\text{me}}{\mathbf{E}}_{\text{avg}}(\omega )+{\overline{\overline{\chi }}}_{\text{mm}}{\mathbf{H}}_{\text{avg}}(\omega )$$

$$\hat{\mathbf{n}}\times \Delta \mathbf{H}=j\omega ϵ{\overline{\overline{\chi }}}_{\text{ee}}{\mathbf{E}}_{\text{avg}}+j\omega \sqrt{\mu ϵ}{\overline{\overline{\chi }}}_{\text{em}}{\mathbf{H}}_{\text{avg}}$$

$$\hat{\mathbf{n}}\times \Delta \mathbf{E}=-j\omega \sqrt{\mu ϵ}{\overline{\overline{\chi }}}_{\text{me}}{\mathbf{E}}_{\text{avg}}-j\omega \mu {\overline{\overline{\chi }}}_{\text{mm}}{\mathbf{H}}_{\text{avg}}.$$

$$\mathbf{E}(\mathbf{r})=$$

$$-\nabla \times {\iint }_{S}G(\mathbf{r},{\mathbf{r}}^{\prime }){\mathbf{K}}^s({\mathbf{r}}^{\prime })d{\mathbf{r}}^{\prime }$$

$$\mathbf{H}(\mathbf{r})=$$

$$+\nabla \times {\iint }_{S}G(\mathbf{r},{\mathbf{r}}^{\prime }){\mathbf{J}}^s({\mathbf{r}}^{\prime })d{r}^{\prime }$$

$$\mathbf{E}(\mathbf{r})$$

$$\mathbf{H}(\mathbf{r})$$

$$(\mathcal{L}\mathbf{X})(\mathbf{r})$$

$$(\mathcal{R}\mathbf{X})(\mathbf{r})$$

$${\nabla }^{2}G(\mathbf{r},{\mathbf{r}}^{\prime })+k^{2}G(\mathbf{r},{\mathbf{r}}^{\prime })=-\delta (\mathbf{r},{\mathbf{r}}^{\prime })$$

$$G(\mathbf{r},{\mathbf{r}}^{\prime })$$

$$G(\mathbf{\rho },{\mathbf{\rho }}^{\prime })$$

$$(\mathcal{L}\mathbf{X})({\mathbf{r}}_i)$$

$$=\sum _{j=1}^N{S}_{i,j}^{\mathcal{L}}\mathbf{X}({\mathbf{r}}_j)$$

$$(\mathcal{R}\mathbf{X})({\mathbf{r}}_i)$$

$$=\sum _{j=1}^N{S}_{i,j}^{\mathcal{R}}\mathbf{X}({\mathbf{r}}_j),$$

$$S_{i,j}^{\mathcal{L}}$$

$$S_{i,j}^{\mathcal{R}}$$

$${\mathbb{E}}^s=-j\omega \mu {\mathbb{S}}^{\mathcal{L}}{\mathbb{J}}^s-{\mathbb{S}}^{\mathcal{R}}{\mathbb{K}}^s$$

$${\mathbb{H}}^s=-j\omega ϵ{\mathbb{S}}^{\mathcal{L}}{\mathbb{K}}^s+{\mathbb{S}}^{\mathcal{R}}{\mathbb{J}}^s$$

$$-\hat{\mathbb{N}}\left({\mathbb{E}}_{s,2}-{\mathbb{E}}_{s,1}\right)$$

$$\hat{\mathbb{N}}\left({\mathbb{H}}_{s,2}-{\mathbb{H}}_{s,2}\right)$$

$${\mathbb{E}}_{\text{so}}={\mathbb{E}}_0,{\mathbb{H}}_{\text{so}}={\mathbb{H}}_0$$

$$\hat{\mathbb{N}}{\mathbb{E}}^s$$

$$\hat{\mathbb{N}}{\mathbb{H}}^s$$

$$\displaystyle\left[\begin{array}[]{cccc}&&\mathbb{I}&\varnothing\\ \lx@intercol\hfil\smash{\raisebox{6.0pt}{\large$\left[j\omega\mu\mathbb{S}^{\mathbf{\mathcal{L}}}\right]$}}\hfil\lx@intercol&0&\mathbb{I}\\ \varnothing&\varnothing&\mathbb{\hat{N}}&\varnothing\\ \varnothing&\varnothing&\varnothing&\mathbb{I}\\ \end{array}\right]\left[\begin{array}[]{c}\mathbb{J}^{s}_{\text{so}}\\ \mathbb{J}^{s}_{\text{pec}}\\ \mathbb{E}^{s}_{\text{pec}}\\ \mathbb{E}^{s}_{\text{so}}\end{array}\right]=\left[\begin{array}[]{c}\varnothing\\ \varnothing\\ \varnothing\\ \mathbb{E}_{0}\end{array}\right]$$

$$\mathbb{SY}=\mathbb{B},$$

$${\hat{\mathbf{n}}}_i\times ({\mathbf{H}}_{s,1}^i-{\mathbf{H}}_{s,2}^i)-\frac{j\omega ϵ{\chi }_{\text{ee,i}}}{2}\left({\mathbf{E}}_{s,1}^i+{\mathbf{E}}_{s,2}^i\right)=0$$

$${\hat{\mathbf{n}}}_i\times ({\mathbf{E}}_{s,1}^i-{\mathbf{E}}_{s,2}^i)+\frac{j\omega {\mu }_0{\chi }_{\text{mm,i}}}{2}\left({\mathbf{H}}_{s,1}^i+{\mathbf{H}}_{s,2}^i\right)=0,$$

$$\hat{\mathbb{N}}({\mathbb{H}}_{\text{m}}^1-{\mathbb{H}}_{\text{m}}^2)-{\mathbb{X}}_e\left({\mathbb{E}}_{\text{m}}^1+{\mathbb{E}}_{\text{m}}^2\right)=0$$

$$\hat{\mathbb{N}}({\mathbb{E}}_{\text{m}}^1-{\mathbb{E}}_{\text{m}}^2)+{\mathbb{X}}_m\left({\mathbb{H}}_{\text{m}}^1+{\mathbb{H}}_{\text{m}}^2\right)=0.$$

Full text

I Introduction

Metasurfaces are 2D counterparts of more general volumetric metamaterials[1]. They are composed of 2D arrays of sub-wavelengths unit cells whose microscopic response can be tailored to engineer the macroscopic response of the metasurfaces. By controlling these susceptibilities, various sophisticated wave transformations can be achieved. With recent intense activity in this area wide variety of metasurfaces have been developed enabling versatile wave transforming applications ranging from Radio Frequencies (RF) to optics [2, 3, 4, 5, 6].

An important problem in metasurface research is the development of fast, efficient and reliable simulation platforms. While the metasurfaces are composed of sub-wavelength resonating cells, they themselves are typically large compared to the wavelengths of operation, so that their numerical simulation is essentially a multi-scale problem. To address this issue, the representation of a physical metasurface is transformed to an ideal zero-thickness model, which is expressed using tensorial effective surface susceptibility densities, $\bar{\bar{\chi}}$ to account for various Electromagnetic (EM) effects, including bi-anisotropy [7]. An equivalent zero thickness model of the metasurface represents a spatial discontinuity and thus is treated using Generalized Sheet Transition Conditions (GSTCs) [8, 9, 10]. Based on GSTCs and the surface susceptibilities, various numerical techniques have been proposed recently to solve for the scattered fields from the metasurface, for both frequency [11, 12] and time-domain [13, 14, 15] simulations using Finite Difference (FD) methods, Finite Element Methods (FEM) [16] and Integral Equations (IE) in the Spectral Domain (SD) [17].

All these numerical methods have been demonstrated for computing the scattered fields from standalone metasurfaces, and are practically suitable for finite regions of spaces only. For cases where the metasurface is placed with various other scatterers as part of an electrically large system these methods become computationally challenging requiring substantial memory and computational resources. To solve such problems more generally, Boundary Element Methods (BEM) have been developed (producing a vast body of literature) which solve the scattered EM fields in terms of physical and equivalent electric and magnetic surface currents in a given volume of interest using integral form of the Maxwell’s equations [18, 19, 20, 21, 22]. As the BEM method does not require meshing of the entire volume, it is computationally efficient and thus well suited for solving electrically large problems.

In this work, an idealized model of the zero thickness metasurface in conjunction with the GSTCs is treated as a generalized boundary condition connecting different volumetric regions, and integrated into the BEM to solve for total scattered fields in electrically large computational domains. Given the generalized field transformation properties of the metasurface, there are non-zero electric and magnetic surface currents that exist on the metasurface. Compared to conventional BEM methods based on Electric Field Integral Equations (EFIE) and Magnetic Field Integral Equations (MFIE) which typically involve only electric surface currents the proposed method thus solves for the total fields in the presence of both electric and magnetic equivalent currents on the metasurface.

The paper is structured as follows. Sec. II outlines the general problem consisting of arbitrary number of finite space regions connected through various conventional boundaries in addition to a set of metasurface boundaries. Sec. III presents the BEM procedure for discretizing Maxwell’s equations. Sec. IV integrates the metasurface boundaries into the BEM, and provides an illustrative example of two regions connected by a single metasurface. Various numerical results are presented in Sec. V followed by conclusions in Sec. VI.

II Problem Formulation

II-A General Scattering Problem & Conventional Boundary Conditions

Consider a general field scattering problem illustrated in Fig. 1, consisting of several volumetric regions of space. Each of the regions of interest are connected to their neighboring regions through various EM boundary conditions, which could represent either physical or purely mathematical boundaries in space. A known source is next applied at one (or more) boundaries in terms of electric and magnetic surface currents, $\mathbf{J_{0}},~{}\mathbf{K_{0}}$ for instance, which produces scattered fields throughout all the regions. The objective here is to compute the total scattered fields in various regions, satisfying all the boundary conditions and solving the Maxwell’s equations self-consistently.

Two boundaries commonly encountered are Perfect Electric Conductors (PECs) and Perfect Magnetic Conductors (PMCs), which are impenetrable boundaries where the total EM fields goes to zero. For instance, in case of a PEC boundary, the tangential E-fields are continuous and zero, while the H-fields are discontinuous giving rise to surface electric currents $\mathbf{J}_{s}$. Similarly for a PMC boundary, the tangential H-fields are continuous and zero, while the E-fields are discontinuous giving rise to surface magnetic currents $\mathbf{K}_{s}$. Formally,

[TABLE]

for PEC and PMC boundaries, respectively, where $\mathbf{H}_{s,1}$ and $\mathbf{E}_{s,1}$ are the total H- and E-fields on the boundary in the incident region of the boundaries and $\mathbf{\hat{n}}_{i}$ is the unit normal vector to the surface at given point.

For other general boundaries for which the tangential E- and H-fields are both discontinuous (i.e. penetrable boundaries), both electric and magnetic surface currents exist and are given by

[TABLE]

where a special case is the interface between two dielectric materials where both $\mathbf{K}_{s}$ and $\mathbf{J}_{s}$ are zero. Completely general boundaries with non-zero $\mathbf{K}_{s}$ and $\mathbf{J}_{s}$ can be realized using Electromagnetic Metasurfaces and described in terms of effective surface polarization densities as described next.

II-B Metasurface Boundaries - Generalized Sheet Transition Conditions (GSTCs)

A metasurface is a two dimensional array of sub-wavelength electromagnetic scatterers with zero thickness ($\delta=0$) which produces a spatial discontinuity in the amplitude and phase of an incoming electromagnetic wave, as shown in Fig. 2. The Generalized Sheet Transition Conditions (GSTCs) were developed by Idemen in [9] to model such these discontinuities and were later applied to metasurfaces in [8]. For a general metasurface embedded inside a uniform media with $(\epsilon,\mu)$, the GSTCs can be written in the frequency-domain as:

[TABLE]

where $\Delta\psi=\psi_{2}-\psi_{1}$ represents the difference between the fields across the metasurface, and $\mathbf{P}$ and $\mathbf{M}$ are the electric and magnetic surface polarization densities. The term $X_{||}$ is the component that is tangential to the metasurface and the term $X_{n}$ is perpendicular to the metasurface. The surface polarization densities are produced in response to a field interacting with the metasurface. These polarizations can be related to the average electric and magnetic fields through the use of surface susceptibilities, and are expressed in general as:

[TABLE]

where $\mathbf{E}_{\text{avg}}=(\mathbf{E}_{s,2}+\mathbf{E}_{s,1})/2$ and $\mathbf{H}_{\text{avg}}=(\mathbf{H}_{s,2}+\mathbf{H}_{s,1})/2$ are the average tangential electric and magnetic field across the metasurface respectively, expressed in terms of the total fields in each of the two regions across the boundaries. $\overline{\overline{\chi}}_{\text{ee}}$ and $\overline{\overline{\chi}}_{\text{mm}}$ are the effective electric and magnetic surface susceptibilities respectively, and $\overline{\overline{\chi}}_{\text{em}}$ and $\overline{\overline{\chi}}_{\text{me}}$ are the cross-anisotropic surface susceptibilities of the metasurface. Let us assume for simplicity that $\mathbf{P}_{n}=\mathbf{M}_{n}=0$ which simplifies Eq. (3) and yields:

[TABLE]

The surface susceptibilities thus set the relationship between all the scattered fields across the metasurface, which can alternatively be synthesized to transform specified incident fields into desired transmission and reflection fields, i.e. total scattered fields [23]. Therefore, the surface susceptibility description of metasurfaces represents a powerful platform to describe arbitrary boundary conditions.

II-C Scattering Formulation

All the boundary conditions above relate the surface currents to the fields just across the boundaries. However, the general goal of the scattering problem illustrated in Fig. 1 is determining the total scattered fields anywhere in the entire computational region. The approach taken here is to determine the fields within each region using an integral representation of Maxwell’s equations in the frequency domain. Within each region the fields are a consequence of the surface currents ($\mathbf{K}_{s}$ and $\mathbf{J}_{s}$) present on that regions boundaries. For impenetrable boundaries the surface currents are determined by a boundary condition such as a PEC or PMC. For penetrable boundaries (interfaces) the currents are such that the interface conditions are maintained. As interfaces allow coupling between regions, the surface currents on these boundaries contribute to the fields in both regions. Some regions will have external boundaries that extend to infinity and allow for free radiation. These boundaries are handled naturally by the integral representation of the electromagnetic equations due to the use of Green’s function that goes to zero at infinity.

Below we will show how, using the BEM method, the integral equations in each region are coupled through the interfaces to form a complete set of self-consistent linear equations that can be solved for the surface currents present for all regions. It should be noted that the resulting surface currents are not typically physical but mathematical artifacts that enforce the boundary conditions and capture the geometrical implications of the field configuration for each region.

III Boundary Element Method (BEM)

III-A Integral Equations for the Regions

When applied to EM the Boundary Element Method (BEM) uses an integral representation of Maxwell’s equations to determine the scattered fields inside of a region. It is assumed that electromagnetic fields are produced by electric and magnetic surface current densities, $\mathbf{J}$ and $\mathbf{K}$, present on the surfaces enclosing a uniform volume of space. These surface current densities are integrated over the entire surface using the frequency domain version of Maxwell’s equations, giving [24, 25]:

[TABLE]

where $\mathbf{E}(\mathbf{r})$ and $\mathbf{H}(\mathbf{r})$ are the electric and magnetic fields inside of the region enclosed by the integration, $k$ is the wave vector in the region, $G$ is the Green’s function for electrodynamics and $\mathbf{J}^{s}$ and $\mathbf{K}^{s}$ are the electric and magnetic surface current densities on the surface, and $\mathbf{r}=(x,y,z)$. Time convention used here is $e^{j\omega t}$. Primed and unprimed variables refer to source and observation locations, respectively. Because Eq. (6a) and (6b) have similar terms, it is useful to represent these equations using a sum of linear operators acting on the surface currents as

[TABLE]

where the operators $\mathbf{\mathcal{L}}$ and $\mathbf{\mathcal{R}}$ are written as:

[TABLE]

These operators involve the use of a Green’s function which characterizes the impulse response of an inhomogeneous linear differential equation [26]. In the case of electrodynamics, the Green’s function is a solution to the Helmholtz equation:

[TABLE]

where $\delta$ is the delta function in space. The use of the Green’s function in Eq. (7) and (8) represents the generation of the fields from the surface currents present at the interfaces. The linearity of the Maxwell’s equations is next used to form the total field response from a superposition of various impulse responses where the source of each impulse is prescribed by the surface current distributions. The solution to Eq. (9) is well known and in two and three dimensions can be written as [24]:

[TABLE]

where $\boldsymbol{\rho}=(x,y)$ is the position vector in two dimensions, and $H_{0}^{(2)}$ is the Hankel function of the first and second kind. Note that the natural boundary condition for both of these functions is that the field is zero at infinity.

For situations with multiple connected regions as shown in Fig. 1, the integral equation operators (Eq. 8) can be applied to each region and interface equations are used to couple the regions. Each region is “extracted” and the scattered fields due to the corresponding surface currents are next computed by operators $\mathbf{\mathcal{L}}^{i}$ and $\mathbf{\mathcal{R}}^{i}$ (for the i’th regions). These surface currents present on the interfaces will be unknowns determined by the complete self-consistent solution of the entire domain.

It should be noted that within each domain the fields are created by currents present on the boundaries of that region only. These currents represent the fields coupled in from the adjacent regions as determined by the interface equations or the imposition of a boundary condition. As such they are fictitious currents that: 1) enforce the interface/boundary conditions; 2) capture the geometry of the region; and 3) take into account the influence of the fields of the surrounding region.

III-B Discretization – BEM with Pulse Functions

The BEM is a well known and thoroughly researched method that calculates the scattered fields interacting with surfaces prescribing a region by descritizing Eq. (7) [24, 25]. To model a region we assume that the surfaces of the region are discritized into a collection of elements – where for our 2D case we would describe the surfaces/interfaces by line segments. In this paper it is assumed that there are two induced electric and magnetic surface currents ($\mathbf{J}_{s}$ and $\mathbf{K}_{s}$) present in the element; both of which “flow” parallel to the surface and are pulse functions and uniform over the entire element[24].111There are implementations of the BEM that use higher order interpolation for the surface currents over the boundary to improve the accuracy of the method [25, 18], however, for simplicity we use uniform elements.

The surfaces prescribing the region are described by a set of $N$ line segments each centered at $\mathbf{r}_{j}$, with a surface normal $\mathbf{\hat{n}}_{j}$ and a length $d\mathbf{r}_{j}$ for the $j$’th segment. Under these assumptions Eq. (8a) and (8b) are descritized by evaluating the Green’s function at the center of the line segment ($\mathbf{r}_{j}$) and using the length $d\mathbf{r}_{j}$ as the weight of the contribution to the sum giving:

[TABLE]

with

[TABLE]

These equations can then be used with Eq. (7) to relate the set of surface fields present at the elements ($\mathbb{E}^{s}$ and $\mathbb{H}^{s}$) to the set of surface currents ($\mathbb{J}^{s}$ and $\mathbb{K}^{s}$) by two matrix equations,

[TABLE]

where $\mathbb{S}^{\mathbf{\mathcal{L}}}$ and $\mathbb{S}^{\mathbf{\mathcal{R}}}$ are matrices that are formed using the $S^{\mathbf{\mathcal{L}}}_{i,j}$ and $S^{\mathbf{\mathcal{R}}}_{i,j}$ coefficients defined in Eq. (12). It is evident that as we are determining the surface fields from the surface currents there is a contribution to each field from the self-same element. As the Green’s function has a singularity at the source location when $\mathbf{r}=\mathbf{r}^{\prime}$ this contribution needs to be handled carefully. Standard procedures exist for extracting this singularity and evaluating this contribution to the total field solution[24].

These equations thus relate the fields created to a known set of surface currents, however, in general (except for a defined source) the currents at the interfaces are unknown and constrained by various interface boundary conditions (BC) discussed in Sec. II. The relationship between the scattered fields and the surface currents represented by Eq. (13) is an additional requirement on the solution that the fields are a solution to Maxwell’s equations. The BEM method thus combines the various BC equations at every interface with Eq. (13) to solve for the unknown surface current distributions, from which the general scattered fields can finally be calculated.

To illustrate this method, let us consider an arbitrarily shaped boundary described in terms of the electric and magnetic surface currents. These boundary conditions can be expressed in matrix form, once the interface is discretized so that

[TABLE]

where $\mathbb{\hat{N}}$ is a matrix operator formed from the operator $(\mathbf{\hat{n}}_{i}\times)$. In addition to these BC’s an analysis will need to define a source. This can be done by defining a sub-set of the interface elements to be a source and prescribing a known electric and magnetic field distribution (or alternatively the electric and magnetic surface currents) on this portion of the surface,

[TABLE]

For many EM problems it is not needed to solve simultaneously for the $\mathbf{E}$ and $\mathbf{H}$ fields. Such a formulation is known as Electric Field Integral Equation (EFIE) which calculates the radiated electric fields obtained from the induced surface currents [24] [27]. In such cases, the EFIE would enforce Eq. (14) while assuming that $\mathbb{K}_{s}=0$. Both Eq. (14) and Eq. (13a) are then solved self-consistently for the unknown local current $\mathbb{J}^{s}$ created by known incident fields. For example, simple scattering of an EM wave from PEC boundary, as shown in Fig. 3, will not generate any $\mathbf{K}_{s}$ terms and only Eq. 6a needs to be solved in conjunction with the appropriate BCs of the PEC,

[TABLE]

In this case, we can form a set of equations from geometry consisting of a single source with surface currents $\mathbb{J}^{s}_{\text{so}}$, surface fields $\mathbb{E}^{s}_{\text{so}}$ and a single PEC surface with surface currents $\mathbb{J}^{s}_{\text{pec}}$ and fields $\mathbb{E}^{s}_{\text{pec}}$. The complete set of surface currents would be $\mathbb{J}^{s}=[\mathbb{J}^{s}_{\text{so}}\;\mathbb{J}^{s}_{\text{pec}}]$ and the surface fields $\mathbb{E}^{s}=[\mathbb{E}^{s}_{\text{pec}}\;\mathbb{E}^{s}_{\text{so}}]$. Using Eq. (13), (15) and (16), we can form the field matrix equation:

[TABLE]

or

[TABLE]

where $\varnothing$ is a Null vector or matrix and $\mathbb{I}$ is an identity matrix. This is a complete set of linear equations that can be solved for $\mathbb{J}^{s}$ and $\mathbb{E}^{s}$, embedded in $\mathbb{Y}$, from which the fields in the entire region can be calculated using Eq. (7) and the operators given in Eq. (11).

An alternative formulation, the Magnetic Field Integral Equation (MFIE), calculates the radiated magnetic fields from the induced surface currents [24][28]. The MFIE is a Fredholm integral equation of the second kind, compared to the EFIE which is of the first kind, however it is generally limited to closed structures [28]. This formulation would solve the dielectric boundary condition by solving both Eq. (14) and (13b) for the unknown local current $\mathbf{J}_{s}$, once again imposing that $\mathbf{K}_{s}=0$. The two methods can be used using a linear combination to eliminate singularities in the formulation which is known as the Combined Field Integral Equation (CFIE)[18].

It should be noted that for a general metasurface boundary condition, relating arbitrary fields across it and described in terms of tensorial surface susceptibilities, both surface currents $\mathbf{J}_{s}$ and $\mathbf{K}_{s}$ exist on the interface, so that both EFIE/MFIE and CFIE methods cannot be straightforwardly applied to this boundary. Consequently we will develop an appropriate formulation of the metasurface boundary next combining GSTCs and the conventional BEM technique.

IV Metasurface Integration in BEM

IV-A Discretized GSTCs

In order to simulate the behavior of a general metasurface acting as a boundary between two adjoining regions of space, GSTCs of Eq. 3 have to be implemented into the BEM equations. Similar to other boundary conditions, the GSTCs relate the electric and magnetic fields on either side of the boundary with each other through the surface susceptibilities. For the sake of simplicity, let us consider a mono-isotropic surface, so that $\overline{\overline{\chi}}_{\text{em}}=\overline{\overline{\chi}}_{\text{me}}=0$ and purely scalar susceptibilities, i.e. $\overline{\overline{\chi}}=\chi$. For such a surface, the metasurface interface condition of Eq. (5) for the $i$’th element when discretized are,

[TABLE]

where the all the fields involved are the tangential fields only to the surface, which can further be expressed in matrix form as:

[TABLE]

where the subscript, $m$, is used to indicate the terms related to metasurface. These equations assume that the fields exit in all three cardinal directions on the surface of the boundary. However, if the surface element’s normal vector is in the same direction as one of the spatial directions the GSTCs will only provide two valid equations instead, related to the tangential fields only. In these cases, the relevant equation is removed and another equation is introduced which enforces that the surface current density perpendicular to the surface is zero.

IV-B Two Region Implementation with Single Metasurface

To illustrate the introduction of the metasurface into a BEM method, we will present for simplicity the mathematical formulation for two regions only. The first with a single source and the two regions simply connected by a metasurface as shown in Fig. 4. In the first region, the surface currents consist of $\mathbb{J}^{s}_{1}=[\mathbb{J}_{\text{so}},\;\mathbb{J}_{\text{m}}]$ and $\mathbb{K}^{s}_{1}=[\mathbb{K}_{\text{so}},\;\mathbb{K}_{\text{m}}]$ and for the second region where only the metasurface is present we have $\mathbb{J}_{s}^{2}=\mathbb{J}_{\text{m}}$ and $\mathbb{K}_{s}^{2}=\mathbb{K}_{\text{m}}$. For each of the regions, the operators $\mathbb{S}^{\mathbf{\mathcal{L}}}$ and $\mathbb{S}^{\mathbf{\mathcal{R}}}$ are formed noting that $\mathbb{J}_{\text{m}}$ and $\mathbb{K}_{\text{m}}$ are present for both, leading to

[TABLE]

Note that the “open” radiating surface is not included in the discretized surface model as it is at infinity. The GSTC expressed in Eq. (17) are then used to related the field $\mathbb{E}^{1}_{\text{m}}$, $\mathbb{E}^{2}_{\text{m}}$, $\mathbb{H}^{1}_{\text{m}}$ and $\mathbb{H}^{2}_{\text{m}}$, and the source equation (15) to define the source fields $\mathbb{E}_{\text{so}}=\mathbb{E}_{0}$ and $\mathbb{H}_{\text{so}}=\mathbb{H}_{0}$. Placing all these equations in a matrix formulation produces Eq. (18).

Which can now finally be solved for all the unknown surface currents and fields and Eq. (7) can be subsequently used to calculate fields anywhere within the two regions.

V Numerical Demonstration

V-A Simulation Setup

Fig. 4 shows the numerical setup consisting of a metasurface of length $\ell$ located at $x=0$. In order to simplify the simulation, a 2D problem is considered where the field varies only in the $x-y$ plane. An input source consisting of both $\mathbf{E}_{s}$ and $\mathbf{H}_{s}$ fields is applied at $x=-d_{s}$ using a source surface of length $\ell$. This source is configured to create a TE field with a diffracting Gaussian-like profile with the waist at $d_{s}$ and width of $\sigma_{y}$.

It should be noted that the separation between the source and the metasurface, $d_{s}$ does not affect the memory requirements of the simulation since the surfaces are linked together using the scattered field equations Eq. (13). In addition, the source and surface discretization can also be different. This implies that the computation region can be arbitrarily large.

The frequency dependence of metasurface susceptibilities are assumed to follow a Lorentzian distribution, given by

[TABLE]

where $\omega_{p}$, $\omega_{0}$ and $\alpha$ are the plasma frequency, resonant frequency and the loss-factor of the oscillator, respectively, and subscripts $e$ and $m$ denote electric and magnetic quantities. For the chosen operation frequency of 60 GHz, the metasurface size is fixed to $\ell=0.1~{}$m and the source surface’s length is set to $\ell_{s}=0.08~{}$m with a separation of $d_{s}=0.05~{}$m. The metasurface length was chosen to be sufficiently large that the source field at the surface was contained within it.

In our numerical simulation, the discretization of the metasurface and the source was set to $n_{\lambda}=40$ divisions per wavelength. The fields in a rectangular region surrounding the metasurface with dimensions $x=[-d_{s},~{}d_{s}]$ and $y=[-\ell/2,~{}\ell/2]$ and discretization $\Delta x=\lambda/10$ and $\Delta y=\lambda/n_{\lambda}$ were calculated using Eq. (13) after the initial simulation.

One advantage of the BEM method is the ability to create visualizations of the field distributions of either the total field (excitation plus scattered fields) or to visualize these fields independently. For example in the first region, one can calculate the incident field by simply applying Eq. 13 to the currents present in the source surface. Conversely, for either region, Eq. 13 can be applied to the currents in the metasurface and the scattered fields (reflected and transmitted) can be determined. Of course, if in the first region both source and metasurface currents are used, then total fields will be calculated.

V-B Simulation Results

To demonstrate the method, two cases of metasurfaces will be considered: 1) A uniform metasurface with $\chi_{\text{ee}}\neq\chi_{\text{mm}}$, and 2) a non-uniform metasurface with spatially varying $\chi_{\text{ee}}(y)=\chi_{\text{mm}}(y)$.

The uniform metasurface is chosen as the first example, since its scattered fields can be readily determined analytically. For a uniform metasurface excited with an arbitrary shaped input beam $\mathbf{E}_{0}(y)$, at a fixed angular frequency, $\omega$, the transmitted and reflected fields in the spatial frequency $k_{y}$, are given by, [29]

[TABLE]

where $\tilde{\mathbf{E}}(k_{y})$ represents the spatial Fourier transform of $\mathbf{E}(y)$, $\mathbf{E}_{t}$ and $\mathbf{E}_{t}$ are the scattered fields in transmission and reflection. This physically represents the transmitted and reflected field response of the metasurface for a specific $k_{y}$, which corresponds to a specific plane-wave excitation (in the propagation regime). The spatial scattered fields $\tilde{\mathbf{E}}_{t}(y)$ and $\tilde{\mathbf{E}}_{r}(y)$ are obtained using inverse Fourier transforms of the above fields: $\mathbf{E}_{t}(x=0+,y)=\mathcal{F}_{y}^{-1}\left\{\tilde{\mathbf{E}}_{t}(k_{y})\right\}$ and $\mathbf{E}_{r}(x=0_{-},y)=\mathcal{F}_{y}^{-1}\left\{\tilde{\mathbf{E}}_{r}(k_{y})\right\}$. This method can be used to validate the numerical results of the proposed method.

Fig. 5 shows the total scattered fields and the scattered fields just before and after the metasurface. The convergence plots in Fig. 6 show the effect of meshing density on the computed fields, where both the transmitted and reflected fields are clearly seen to be converging to analytical values, which indicates that with higher discretization they should approach even closer to the expected results. This provides a good validation of the method.

Next, a non-uniform metasurface is considered whose electric and magnetic surface susceptibilities are assumed equal but modulated in space, emulating a diffraction grating. For a physical metasurface, this can be achieved by periodically modulating the resonant frequency $\omega_{0}$ of the Lorentzian function of Eq. (20), given by

[TABLE]

where $\Delta_{m}$ controls the intensity of the modulation and $\beta_{m}$ is the spatial frequency of modulation. To produce strong diffraction orders, $\Delta_{m}=0.1$ and $\beta_{m}=k/5$ were used with a Gaussian beam of width $\sigma_{y}=4/\beta_{m}$. Fig. 7a shows the calculated total power in the two regions where the normally-incident Gaussian beam is split into several diverging higher-order diffraction order beams. To better visualize these diffraction orders, the transmitted field of the metasurface is captured and a spatial Fast Fourier Transform (FFT) is applied. Fig. 7b shows the strengths of various diffraction orders as a function of the mesh density, $n_{\lambda}$. Resulting diffractions orders are first seen to be equally spaced with $k_{y}=\beta_{m}$ as expected, and their strengths are gradually converged to a constant value beyond which higher meshing has no impact on the results. To validate these results, the strengths of the generated harmonics are compared with a Yee-cell based Finite Difference Frequency Domain (FDFD) method of [12], using the same metasurface parameters with a discretization of $n_{\lambda}=40$. Fig. 7b also shows this comparison, where an excellent agreement is observed between the FDFD and the proposed method, so that the validity of the method is clearly established.

VI Conclusion

A methodology for determining the scattered EM fields present for interconnected regions with metasurface boundaries has been proposed. The method is based on a BEM formulation of the frequency domain version of Maxwell’s equations – which expresses the fields present in a region due to surface currents on the boundaries. Multiple regions are coupled together by shared surface currents which can be solved for in a self-consistent manner. The general metasurface boundaries represented as surface susceptibilities were integrated in the BEM method next using GSTCs.

To illustrate the method, two examples were presented for EM scattering of a Gaussian beam. Firstly a uniform metasurface, and secondly a spatially modulated metasurface. For the uniform surface the BEM results were compared to a semi-analytical method and it was shown that as the boundary segment length was decreased, the BEM results converged to the semi-analytical results. The second example was for a spatially modulated metasurface (essentially a grating) and the BEM results were compared to a Yee-cell frequency domain method. The two methods were found to predict essentially the same diffraction components confirming the BEM accuracy. It can be noted that the BEM method as formulated used a very simple uniform element function and more sophisticated methods would allow for larger elements to be used. However, for this paper the purpose is to simply show the functionality of the method.

Although, this paper shows a simple application of the BEM method with an incorporated metasurface, it establishes the applicability of the method for larger scale scattering problems. Moreover, the GSTC model of metasurface boundaries may also be seen as pure numerical tool to connect multiple regions by mimicking and transcending the functionalities of conventional boundaries. Situations where the metasurface is a part of a larger problem with multiple scattering objects (both electrically large and small), and scattering response is the prime objective (similar to Radar Cross-Section, RCS), the proposed approach may prove to be an invaluable tool. Such problems have been a traditional domain for BEM modeling and with the appropriate incorporation of advanced BEM methods and computational techniques, the proposed methodology is likely to be increasingly useful.