Generalized Sheet Transition Condition FDTD Simulation of Metasurface

TL;DR

This paper introduces a novel FDTD simulation method based on GSTCs for modeling complex, nonlinear, and dynamic metasurfaces, extending standard FDTD to handle zero-thickness, polychromatic, and space-time varying cases.

Contribution

It presents a new FDTD scheme incorporating GSTCs for metasurfaces, enabling accurate simulation of complex metasurface behaviors with a simple, fully numerical approach.

Findings

Validated with five representative examples.

Can simulate polychromatic and nonlinear metasurfaces.

Easily extendable to bianisotropic cases.

Abstract

We propose an FDTD scheme based on Generalized Sheet Transition Conditions (GSTCs) for the simulation of polychromatic, nonlinear and space-time varying metasurfaces. This scheme consists in placing the metasurface at virtual nodal plane introduced between regular nodes of the staggered Yee grid and inserting fields determined by GSTCs in this plane in the standard FDTD algorithm. The resulting update equations are an elegant generalization of the standard FDTD equations. Indeed, in the limiting case of a null surface susceptibility ($\chi_\text{surf}=0$), they reduce to the latter, while in the next limiting case of a time-invariant metasurface $[\chi_\text{surf}\neq\chi_\text{surf}(t)]$, they split in two terms, one corresponding to the standard equations for a one-cell ($\Delta x$) thick slab with volume susceptibility ($\chi$), corresponding to a diluted approximation ($\chi=\chi_\text{surf}/(2\Delta x)$) of the zero-thickness target metasurface, and the other transforming this slab in a real (zero-thickness) metasurface. The proposed scheme is fully numerical and very easy to implement. Although it is explicitly derived for a monoisotropic metasurface, it may be straightforwardly extended to the bianisotropic case. Except for some particular case, it is not applicable to dispersive metasurfaces, for which an efficient Auxiliary Different Equation (ADE) extension of the scheme is currently being developed by the authors. The scheme is validated and illustrated by five representative examples.

Equations

$$\displaystyle\left(\begin{array}[]{c}-\Delta\tilde{H}_{z}\\ \Delta\tilde{H}_{y}\\ \end{array}\right)$$

$$\displaystyle+j\omega\sqrt{\varepsilon_{0}\mu_{0}}\left(\begin{array}[]{cc}\tilde{\chi}_{\textrm{em}}^{yy}&\tilde{\chi}_{\textrm{em}}^{yz}\\ \tilde{\chi}_{\textrm{em}}^{zy}&\tilde{\chi}_{\textrm{em}}^{zz}\\ \end{array}\right)\left(\begin{array}[]{c}\tilde{H}_{y,\textrm{av}}\\ \tilde{H}_{z,\textrm{av}}\\ \end{array}\right),$$

$$\displaystyle\left(\begin{array}[]{c}-\Delta\tilde{E}_{y}\\ \Delta\tilde{E}_{z}\\ \end{array}\right)$$

$$\displaystyle+j\omega\sqrt{\varepsilon_{0}\mu_{0}}\left(\begin{array}[]{cc}\tilde{\chi}_{\textrm{me}}^{zz}&\tilde{\chi}_{\textrm{me}}^{zy}\\ \tilde{\chi}_{\textrm{me}}^{yz}&\tilde{\chi}_{\textrm{me}}^{yy}\\ \end{array}\right)\left(\begin{array}[]{c}\tilde{E}_{z,\textrm{av}}\\ \tilde{E}_{y,\textrm{av}}\\ \end{array}\right),$$

$$\Delta H_y$$

$$\Delta E_z$$

$$H_y^{n+\frac{1}{2}}\left(i\right)=H_y^{n-\frac{1}{2}}\left(i\right)+\frac{\Delta t}{{\mu }_0\Delta x}\left[E_z^n\left(i+1\right)-E_z^n\left(i\right)\right],$$

$$E_z^n\left(i\right)=E_z^{n-1}\left(i\right)+\frac{\Delta t}{{\epsilon }_0\Delta x}\left[H_y^{n-\frac{1}{2}}\left(i\right)-H_y^{n-\frac{1}{2}}\left(i-1\right)\right]$$

$$E_z^n\left(i_d+1\right)=$$

$$\frac{\Delta t}{{\epsilon }_0\Delta x}\left[H_y^{n-\frac{1}{2}}\left(i_d+1\right)-H_y^{n-\frac{1}{2}}\left(0^{+}\right)\right],$$

$$H_y^{n-\frac{1}{2}}\left(0^{+}\right)=H_y^{n-\frac{1}{2}}\left(i_d\right)+\frac{{\epsilon }_0}{\Delta t}[{\left({\chi }_{\text{ee}}^{zz}{E}_{z,\text{av}}\right)}^n-{\left({\chi }_{\text{ee}}^{zz}{E}_{z,\text{av}}\right)}^{n-1}],$$

$$E_z^n\left(i_d+1\right)$$

$$\frac{\Delta t}{{\epsilon }_0\Delta x}\left[H_y^{n-\frac{1}{2}}\left(i_d+1\right)-H_y^{n-\frac{1}{2}}\left(i_d\right)\right]-$$

$$\frac{{\chi }_{\text{ee}}^{zz,n}}{2\Delta x}{E}_z^n\left(i_d\right)+\frac{{\chi }_{\text{ee}}^{zz,n-1}}{2\Delta x}{E}_z^{n-1}\left(i_d\right),$$

$$A_{\text{ee}}^{zz,n}=1+\frac{{\chi }_{\text{ee}}^{zz,n}}{2\Delta x}.$$

$$H_y^{n+\frac{1}{2}}\left(i_d\right)=H_y^{n-\frac{1}{2}}\left(i_d\right)+\frac{\Delta t}{{\mu }_0\Delta x}\left[E_z^n\left(0^{-}\right)-E_z^n\left(i_d\right)\right],$$

$$E_z^n\left(0^{-}\right)$$

$$\frac{{\mu }_0}{\Delta t}\left[{\left({\chi }_{\text{mm}}^{yy}{H}_{y,\text{av}}\right)}^{n+\frac{1}{2}}-{\left({\chi }_{\text{mm}}^{yy}{H}_{y,\text{av}}\right)}^{n-\frac{1}{2}}\right],$$

$$H_y^{n+\frac{1}{2}}$$

$$\frac{\Delta t}{{\mu }_0\Delta x}\left[E_z^n\left(i_d+1\right)-E_z^n\left(i_d\right)\right]-$$

$$\frac{{\chi }_{\text{mm}}^{yy,n+\frac{1}{2}}}{2\Delta x}{H}_y^{n+\frac{1}{2}}\left(i_d+1\right)+\frac{{\chi }_{\text{mm}}^{yy,n-\frac{1}{2}}}{2\Delta x}{H}_y^{n-\frac{1}{2}}\left(i_d+1\right).$$

$$A_{\text{mm}}^{yy,n}=1+\frac{{\chi }_{\text{mm}}^{yy,n}}{2\Delta x}.$$

$$H_y^{n+\frac{1}{2}}$$

$$\frac{\Delta t}{{\mu }_0{\mu }_r\Delta x}\left[E_z^n\left(i_d+1\right)-E_z^n\left(i_d\right)\right],$$

$$E_z^n$$

$$\frac{\Delta t}{{\epsilon }_0\Delta x}\left[H_y^{n-\frac{1}{2}}\left(i_d+1,j_d\right)-H_y^{n-\frac{1}{2}}\left(0^{+},j_d\right)\right]+$$

$$\frac{\Delta t}{{\epsilon }_0\Delta y}\left[H_x^{n-\frac{1}{2}}\left(i_d+1,j_d\right)-H_x^{n-\frac{1}{2}}\left(i_d+1,j_d\right)\right],$$

$$H_y^{n-\frac{1}{2}}\left(0^{+},j_d\right)=$$

$$\frac{{\epsilon }_0}{\Delta t}\left[{\left({\chi }_{\text{ee}}^{zz}{E}_{z,\text{av}}\right)}^n-{\left({\chi }_{\text{ee}}^{zz}{E}_{z,\text{av}}\right)}^{n-1}\right],$$

$$E_z^n$$

$$\frac{\Delta t}{{\epsilon }_0\Delta x}\left[H_y^{n-\frac{1}{2}}\left(i_d+1,j_d\right)-H_y^{n-\frac{1}{2}}\left(i_d,j_d\right)\right]+$$

$$\frac{\Delta t}{{\epsilon }_0\Delta y}\left[H_x^{n-\frac{1}{2}}\left(i_d+1,j_d\right)-H_x^{n-\frac{1}{2}}\left(i_d+1,j_d\right)\right]-$$

$$\frac{1}{2\Delta x}\left[{\chi }_{\text{ee}}^{zz,n}{E}_z^n\left(i_d,j_d\right)-{\chi }_{\text{ee}}^{zz,n-1}{E}_z^{n-1}\left(i_d,j_d\right)\right].$$

$$A_{\text{ee}}^{zz,n}=1+\frac{{\chi }_{\text{ee}}^{zz,n}}{2\Delta x},$$

$$\begin{array}{rl}& H_y^{n+\frac{1}{2}}\left(i_d,j_d\right)=H_y^{n-\frac{1}{2}}\left(i_d,j_d\right)+\\ & \frac{\Delta t}{{\mu }_0\Delta x}\left[E_z^n\left(0^{-},j_d\right)-E_z^n\left(i_d,j_d\right)\right],\end{array}$$

$$E_z^n\left(0^{-},j_d\right)$$

$$\frac{{\mu }_0}{\Delta t}\left[{\left({\chi }_{\text{mm}}^{yy}{H}_{y,\text{av}}\right)}^{n+\frac{1}{2}}-{\left({\chi }_{\text{mm}}^{yy}{H}_{y,\text{av}}\right)}^{n-\frac{1}{2}}\right],$$

$$H_y^{n+\frac{1}{2}}$$

$$\frac{\Delta t}{{\mu }_0\Delta x}\left[E_z^n\left(i_d+1,j_d\right)-E_z^n\left(i_d,j_d\right)\right]-$$

$$\frac{{\chi }_{\text{mm}}^{yy,n+\frac{1}{2}}}{2\Delta x}{H}_y^{n+\frac{1}{2}}\left(i_d+1,j_d\right)+$$

$$\frac{{\chi }_{\text{mm}}^{yy,n-\frac{1}{2}}}{2\Delta x}{H}_y^{n-\frac{1}{2}}\left(i_d+1,j_d\right).$$

$$\Delta {\tilde{H}}_y$$

$$\Delta {\tilde{E}}_z$$

Full text

I Introduction

Metasurfaces are two-dimensional subwavelengthly thin bianisotropic structures consisting of scattering particles distributed on a substrate. They can perform a great diversity of electromagnetic wave transformations such as, for instance, broad-band focusing [1], hologram generation [2], anomalous refraction and reflection [3], orbital angular momentum production [4], space-surface wave manipulations [5] and remote processing [6].

Figure 1 represents a general metasurface, transforming an incident wave into reflected and transmitted waves. As indicated in the figure, a metasurface is typically much thinner than the operation wavelength, i.e. $\delta\ll\lambda$. Therefore, it can be modeled as a zero-thickness sheet, which may simultaneously represent an electric and magnetic discontinuity. The metasurface may generally exhibit bianisotropy and variation in both space and time [7]. The only known approach for exactly modeling a general metasurface is the Generalized Sheet Transition Condition (GSTC) technique [8, 9, 10, 11, 6].

Unfortunately, there is currently no commercial software including GSTCs and hence capable to simulate general metasurfaces. A metasurface may be approximated by an electrically thin slab with a volume susceptibility equal to the metasurface susceptibility divided by the thickness of the slab [12, 13]. However, this technique is often inaccurate because a slab is a two-interface structure that cannot reduce to a single sheet structure. For example, as shown in [12, 13], the thin-slab approximation erroneously restricts the absorption level of a perfectly absorbing metasurface to $87\%$.

In [14, 15], a Finite-Difference Time-Domain (FDTD) based technique is proposed for the simulation of single magnetic field and single electric field sheet discontinuities corresponding for instance to graphene and 2D electron gas problems, respectively. In [13] and [16], GSTC-based Finite difference Frequency Domain (FDFD) and Finite Element Method (FEM) schemes are introduced to simulate monochromatic, linear and time-invariant general-discontinuity (electric and magnetic) metasurfaces, respectively. While these approaches generally require meshing the three dimensional space, a GSTC-based spectral-domain integral-equation technique is proposed in [17] to reduce meshing to the surface of the metasurface. An equivalent RLC circuit based FDTD approach is presented in [18, 19] to simulate space-time varying metasurfaces; however, this approach is limited to Lorentzian dispersion and inapplicable to problems involving extra scatterers.

In order to overcome the aforementioned limitations, and hence simulate general polychromatic, nonlinear and space-time varying metasurfaces placed in arbitrary scattering environments, we introduce here a completely general GSTC-FDTD scheme, that is easy to implement and that may be readily integrated in commercial softwares. This scheme is also shown to provides a fundamental generalization of the conventional FDTD technique.

The paper is organized as follows. Section II briefly describes the metasurface synthesis equations in the time domain and the frequency domain. Upon this basis, Sec. III derives the proposed GSTC-FDTD scheme for 1D and 2D problems. Benchmark and illustrative examples are provided in Sec. VI. Finally, conclusions are given in Sec. VII.

II Metasurface Synthesis Equations

The metasurface is assumed to be positioned in the $y-z$ plane of a cartesian coordinate system as shown in Fig. 1. It is illuminated by an arbitrary wave from the left side and transforms this incident wave into arbitrary reflected and transmitted waves at the left and right sides of the structure, respectively. The metasurface is best modelled by electric and magnetic surface susceptibility tensors in the frequency domain or the time domain.

II-A Frequency Domain

Without loss of generality and for simplicity, we assume that the metasurface does not support any normal electric and magnetic polarization currents, so that the corresponding polarization densities are zero, i.e. $P_{x}=0$ and $M_{x}=0$. In this case, the frequency-domain (time harmonic dependency $e^{+j\omega t}$) bianisotropic GSTCs read [20]

[TABLE]

where $\Delta\tilde{\psi}_{u}=\tilde{\psi}_{u}^{\textrm{tr}}-(\tilde{\psi}_{u}^{\textrm{ref}}+\tilde{\psi}_{u}^{\textrm{inc}})$ and $\tilde{\psi}_{u,\textrm{av}}=\frac{(\tilde{\psi}_{u}^{\textrm{inc}}+\tilde{\psi}_{u}^{\textrm{ref}})+\tilde{\psi}_{u}^{\textrm{tr}}}{2}$ with $u=x,y$ or $z$ and $\tilde{\psi}$ representing the spectral electric or magnetic fields, and the superscripts ”inc”, ”tr” and ”ref” denote the incident, transmitted and reflected waves, respectively, and the $\tilde{\chi}$’s represent the generally frequency-dependent surface susceptibility tensor terms.

The synthesis procedure consists in determining the susceptibility tensors in (1) for specified sets (possibly multiple ones for multiple transformations) of fields $(\psi_{u}^{\textrm{inc}},\psi_{u}^{\textrm{ref}},\psi_{u}^{\textrm{tr}})$ [20]. Since Eqs. (1) represent a linear system of equations, this generally (always in the case of a full-rank system) leads to explicit solutions, $\bar{\bar{\chi}}_{\textrm{ab}}=\bar{\bar{\chi}}_{\textrm{ab}}(y,z)$ (a,b$=$e,m). Moreover, if the fields are specified in closed-form, then the solution expressions are also closed form.

II-B Time Domain

The time-domain GSTCs are obtain by taking the inverse Fourier transform of (1). Assuming that the metasurface is non-dispersive, i.e. that the $\tilde{\chi}$’s in (1) do not depend on $\omega$, the terms $\tilde{\chi}\tilde{\psi}_{u}(\omega)$ simply transform to $\chi\psi_{u}(t)$ where $\chi=\tilde{\chi}$, without involving any convolution. The time-domain counterpart of (1) is then obtained by just replacing $j\omega$ with the operator $\frac{d}{dt}$. For notational simplicity and without loss of generality, we write here only the monoisotropic version of the resulting equations:

[TABLE]

Note that $\frac{d}{dt}$ generally operates on the products of the susceptibility tensors and average fields. Such products, when both terms are time-dependent, lead to the generation of new frequencies. Here, we are generally looking for solutions of the form $\bar{\bar{\chi}}_{\textrm{ab}}(y,z;t)$. In contrast to the frequency-domain case, these solutions are generally neither explicit nor closed-form. So, the synthesis problem is substantially more complicated in the time domain than in the frequency domain.

III GSTC-FDTD Scheme

We shall restrict to the case of monoisotropic metasurfaces. We first address the 1D problem, where the metasurface reduces to a point, and is hence a 0D metasurface or “metapoint”, and next extend the scheme to the corresponding 2D problem, where the metasurface reduces to a segment, and is hence a 1D metasurface or “metasegment”. Upon this foundation, the extension to the bianisotropic and 3D cases is straightforward although tedious.

III-A 1D Analysis

We consider here the 1D computational problem of a 0D metasurface, represented in Fig. 2. This still represent the practical 3D problem of a uniform 2D metasurface perpendicularly illuminated by a TEM plane wave. The metasurface is positioned in the 1D FDTD staggered Yee grid as shown in Fig. 2, i.e. at $x=0$, between the two neighboring cells $i_{d}$ and $i_{d}+1$, where $i_{d}$ is the number of the cell before the metasurface. The incident wave illuminates the metasurface from the left of the metasurface ($x<0$) while the right side ($x>0$) is the transmitted-wave region.

Everywhere except at the nodes around the discontinuity, one can use the conventional FDTD update equations. These equations read [21]

[TABLE]

where $n$ represents the discrete time, $t=n\Delta t$, with $\Delta t$ being the time step.

In contrast, special treatment must be applied to the nodes around the discontinuity. Consider the update equation for $E_{z}^{n}\left(i_{d}+1\right)$ in Fig. 2, that would conventionally depend on the surrounding fields $H_{y}^{n-1/2}\left(i_{d}\right)$ and $H_{y}^{n-1/2}\left(i_{d}+1\right)$. Because of the presence of the discontinuity between cells $i_{d}$ and $i_{d}+1$, $H_{y}^{n-1/2}\left(i_{d}\right)$ cannot be readily used as such in (3b) for otherwise nothing would be accounting for the effect of the discontinuity. To remedy this situation, we introduce a virtual magnetic node (H-node) at the right side of the metasurface and rewrite the update equation involving this new node as

[TABLE]

where $H_{y}^{n-\frac{1}{2}}\left(0^{+}\right)$ is unknown. To determine this quantity, we invoke the GSTC equation (2a), whose time-domain discretization reads

[TABLE]

where the effect of the metasurface is now accounted for in the second term of the right-hand side of this equation. Substituting (5) in (4), and grouping identical components, yields

[TABLE]

Now consider the update equation for $H_{y}^{n+\frac{1}{2}}\left(i_{d}\right)$, that would conventionally depend on the surrounding fields $E_{z}^{n}\left(i_{d}\right)$ and $E_{z}^{n}\left(i_{d}+1\right)$. Following the same argument as for $E_{z}^{n}\left(i_{d}+1\right)$, we introduce the virtual electric node (E-node) at the left side of the metasurface and rewrite the update equation (3a) involving this new node as

[TABLE]

next invoke the GSTC equation (2b), whose time-domain discretization reads

[TABLE]

and finally substitute (8) into (7) to obtain

[TABLE]

Equations (6) and (9) are the update equations for the nodes around the metasurface discontinuity. In the absence of metasurface, corresponding to $\chi_{\textrm{ee}}^{zz,n}=\chi_{\textrm{mm}}^{yy,n}=0$, these equations reduce to their conventional counterpart (3), as expected. A next limit case is that of a time-invariant metasurface, for which $\chi_{\textrm{ee}}^{zz,n}$ and $\chi_{\textrm{mm}}^{yy,n}$ are independent of time. In this case, Eq. (9), for instance, reduces to

[TABLE]

with $\mu_{\textrm{r}}=1+\frac{\chi_{\textrm{mm}}^{xx}}{2\Delta x}$. This equation corresponds to its conventional counterpart for a “metasurface” modeled by a thin slab of volume permeability $\mu_{r}$, which may be understood as a dilution of the metasurface surface permeability, $\mu_{\textrm{r,{surf}}}$, across two grid cells, with the overall thickness of $2\Delta x$, namely $\mu_{\textrm{r}}=\mu_{\textrm{r,{surf}}}/(2\Delta x)$. A similar time-invariant susceptibility result is naturally obtained for (6). Thus, the last two terms in (6a) and (9a) may be interpreted as corrections of the corresponding conventional equations to model a zero-thickness sheet.

III-B 2D Analysis

We now consider the 2D computational problem of a 1D metasurface, represented in Fig. 3. This still represents the practical 3D problem of a 1D-nonuniform (e.g. phase-gradient) 2D metasurface illuminated by an oblique TMz or TEz plane wave. Without loss of generality, we solve here the TMz problem, with $E_{z},H_{x}$ and $H_{y}$ being the non-zero field components, while a similar scheme straightforwardly applies to the TEz case. The metasurface is positioned in the 2D FDTD staggered Yee grid as shown in Fig. 3, where $i_{d}$ is the number of the cell before the metasurface and $j_{d}=n_{l}:n_{h}$ represents the metasurface extension in the $y-$direction. The incident wave illuminates the metasurface in the $xz-$ plane.

Consider the update equation for $E_{z}^{n}\left(i_{d}+1,j_{d}\right)$ at the nodes around the discontinuity which would conventionally require $H_{y}^{n-\frac{1}{2}}\left(i_{d},j_{d}\right)$. However, as in the 1D case, this field quantity cannot readily be used as such, because that would miss the effect of the metasurface, and we therefore introduce a virtual array of magnetic nodes at the lower side of the discontinuity. This leads to the update equation

[TABLE]

which is the 2D counterpart of (4). In this expression, $H_{y}^{n-\frac{1}{2}}\left(0^{+},j_{d}\right)$ is determined from the GSTC equation (2a) as

[TABLE]

which is the 2D counterpart of (5). Finally, substituting (12) into (11) and regrouping yields

[TABLE]

the 2D counterpart of (6).

Similarly, one obtains for the 2D counterpart of (7), (8) and (9), assuming a virtual array of electric nodes at the upper side of the discontinuity,

[TABLE]

[TABLE]

and

[TABLE]

IV Note on Dispersive Metasurfaces

Section III, based on (2), assumed non-dispersive time-varying metasurfaces, and the time dependency of the susceptibilities corresponding to the time variance is apparent from the $n$-dependency of the susceptibility terms $\chi_{\textrm{ee}}^{zz,n}$ in (6) and (13) and $\chi_{\textrm{mm}}^{yy,n-1/2}$ in (9) and (16). However, metasurfaces are often dispersive, due to their typically resonant nature, and the corresponding time-varying susceptibilities would then be depending also on frequency, i.e. $\tilde{\chi}=\tilde{\chi}(\omega;t)$ in Sec. II-A. In such a case, the virtual node concept introduced in this work for accommodating GSTCs in FDTD is still applicable, but must be accompanied with a standard FDTD treatment of dispersion, with update equations that are specific to the dispersion considered [22]. A complete extension of the present work to the case of dispersive metasurfaces will be the object of a next paper.

It turns out that the scheme of Sec. III still applies to particular dispersion forms, specifically to those forms that allow proper staggered Yee grid discretization in time. Let us see this in the forthcoming examples.

First consider a metasurface absorbing part of a normally incident wave and transmitting the rest of it. In general, such a metasurface is represented by susceptibilities of the form $\tilde{\chi}_{\textrm{ee}}^{zz}(\omega)=\tilde{\chi}_{\textrm{mm}}^{yy}(\omega)=\frac{\kappa}{j\omega}$, where $\kappa$ only depends on the reflection and transmission coefficients111Assuming TMz operation ($\tilde{H}_{z}=\tilde{E}_{x}=\tilde{E}_{y}=0$) and monoisotropy ($\bar{\bar{\chi}}_{\textrm{em}}=\bar{\bar{\chi}}_{\textrm{me}}=0$, $\bar{\bar{\chi}}_{\textrm{ee}}=\hat{z}\tilde{\chi}_{\textrm{ee}}^{zz}\hat{z}$ and $\bar{\bar{\chi}}_{\textrm{mm}}=\hat{y}\tilde{\chi}_{\textrm{mm}}^{yy}\hat{y}$), the metasurface synthesis equations (1) reduce to $\Delta\tilde{H}_{y}=j\omega\varepsilon_{0}\tilde{\chi}_{\textrm{ee}}^{zz}\tilde{E}_{z,\textrm{av}}$ and $\Delta\tilde{E}_{z}=j\omega\mu_{0}\tilde{\chi}_{\textrm{mm}}^{yy}\tilde{H}_{y,\textrm{av}}$. Inserting the specified fields, $\tilde{E}_{z}^{\textrm{inc}}=e^{-jk_{0}x}$, $\tilde{E}_{z}^{\textrm{ref}}=\Gamma e^{+jk_{0}x}$ and $\tilde{E}_{z}^{\textrm{tr}}=Te^{-jk_{0}x}$, into these equations at $x=0$ yields $\tilde{\chi}_{\textrm{ee}}^{zz}=\frac{1}{j\omega\varepsilon_{0}}\frac{\tilde{H}_{y}^{\textrm{tr}}-(\tilde{H}_{y}^{\textrm{inc}}+\tilde{H}_{y}^{\textrm{ref}})}{[\tilde{E}_{z}^{\textrm{tr}}+(\tilde{E}_{z}^{\textrm{inc}}+\tilde{E}_{z}^{\textrm{ref}})]/2}=\frac{2c_{0}}{j\omega}\frac{T-1-\Gamma}{T+1+\Gamma}$ and $\tilde{\chi}_{\textrm{mm}}^{yy}=\frac{1}{\omega\mu_{0}}\frac{\tilde{E}_{z}^{\textrm{tr}}-(\tilde{E}_{z}^{\textrm{inc}}+\tilde{E}_{z}^{\textrm{ref}})}{[\tilde{H}_{y}^{\textrm{tr}}+(\tilde{H}_{y}^{\textrm{inc}}+\tilde{H}_{y}^{\textrm{ref}})]/2}=\frac{2c_{0}}{j\omega}\frac{T-1-\Gamma}{T+1+\Gamma}$.. Substituting these susceptibilities into (1) yields

[TABLE]

which represents the very particular case where the frequency dependency has vanished. The corresponding discretized time-domain equations read

[TABLE]

From this point, one can follow exactly the same method as in Sec. III with just replacing (5) and (8) for the 1D case and (12) and (15) for the 2D case by (19a) and (19b), respectively.

Consider now the Lorentzian metasurface, $\tilde{\chi}_{ee}^{zz}(\omega)=\frac{1}{-\omega^{2}+j\gamma\omega+\omega_{0}^{2}}$. In this case, without any loss of generality and for the sake of simplicity, we use only $\tilde{\chi}_{ee}^{zz}$, the same result being applicable for the cases involving other components of the susceptibility tensors. Substituting this susceptibility into (1) yields

[TABLE]

and the discretized time-domain counterpart of this equation is obtained by replacing $j\omega$ with $\frac{d}{dt}$, which yields

[TABLE]

This equation requires H-field samples at both half-integer and integer times whereas, according the staggered nature of the Yee grid, the H-field is available only at half-integer times. Equation (21) is thus incompatible with the staggered Yee grid and standard FDTD algorithm. Resolving this issue would require the introduction of auxiliary functions [22] and, as previously mentioned, this will be done in a next paper. Equation (21) would be directly compatible with the FDTD algorithm only in the case $\gamma=0$, which would correspond to a non-causal lossless and hence impractical dispersive system.

Generalizing from the above examples, we see the technique proposed in this paper is directly applicable only when the first and second derivatives of the magnetic or electric field do not simultaneously appear in the time-domain GSTC equation or, equivalently, when $j\omega$ and $(j\omega)^{2}$ do not simultaneously appear in the denominator or nominator of the frequency-domain susceptibility.

V Note on Nonlinear Metasurfaces

The proposed GSTC-FDTD technique is also applicable to nonlinear metasurfaces, where the susceptibility function depends on orders of the electric or magnetic fields larger than one and where the susceptibility tensors reach an order corresponding to the order of the nonlinearity. The nonlinear extension of the proposed technique is straightforward. It has in fact already been implemented in a work by some of the authors of this paper specifically dealing with nonlinear metasurfaces and the details are available in the appendix of [23].

VI Benchmark and Illustrative Examples

This section presents the following five application examples of the GSTC-FDTD scheme presented in Sec. III:

1D analysis of constant-susceptiblity reflection-less 0D metasurface benchmarked with analytical result; 2. 2.

1D analysis of sinusoidally time-varying reflection-less 0D metasurface; 3. 3.

2D analysis of 1D metasurface representing a graphene sheet benchmarked with literature results; 4. 4.

2D analysis of 1D space-time reflection-less varying and dispersive metasurface benchmarked with synthesis specification. 5. 5.

1D analysis of reflection-less half-absorbing and half-transmitting 0D dispersive metasurface benchmarked with synthesis specification;

VI-A Example 1

Consider the 1D computational problem of a 0D reflectionless metasurface with constant and equal susceptibilities $\chi_{\textrm{ee}}^{zz}=\chi_{\textrm{mm}}^{yy}=5$. It may be shown that in the case of equal monoisotropic susceptibilities, the metasurface is reflection-less [20]. The following analytical closed-form transmitted field is found from (2) for the (necessarily) normally incident and transverse electromagnetic (TEM) plane wave $H_{y}^{\textrm{inc}}=\sin\left(\omega t-k_{0}x\right)$:

[TABLE]

The GSTC-FDTD simulation and exact analytical results for this problem are compared in Fig. 4. In the steady state, the two results feature undistinguishable amplitude waveforms in the scale range of the figure, with a zero reflected field, consistently with the design.

VI-B Example 2

The second example is the 1D computational problem of a 0D reflection-less time-varying metasurface with susceptibilities $\chi_{\textrm{ee}}^{zz}=\chi_{\textrm{mm}}^{yy}=\chi=1+\sin\left(\omega t\right)$, where the two susceptibilities have been made equal for zero reflection [20], with normally incident TEM wave $E_{z}=\sin\left(\omega t\right)$. Due to the time-varying nature of the metasurface, generation of new time-harmonic frequencies are expected, as may be understood from inspecting (2).

The simulation result is plotted in Fig. 5. Reflection has been verified to be zero and the emergence of new frequencies may be deduced from the distorted shaped of the transmitted wave in the figure. Unfortunately, this example cannot be benchmarked due to the inexistence of appropriate commercial code.

VI-C Example 3

We shall now demonstrate the applicability of our GSTC-FDTD scheme to a graphene sheet [14] 2D-analyzed as a 1D metasurface. Since the scheme applies to both electric and magnetic field discontinuities, it is expected to automatically work for the case of single magnetic-field discontinuity represented by graphene.

Since graphene is continuous in terms of the electric field (while being discontinuous in terms of the magnetic field), it is naturally positioned at the $E_{z}\left(i_{d}+1,j_{d}\right)$-nodes in the FDTD grid of Fig. 3. We have thus to compute the magnetic field at the nodes before and after the graphene sheet. From [14], the discretized equation for (non-magnetized) graphene is

[TABLE]

with $C^{+}=\frac{\Delta t+2\tau}{2\Delta t\sigma_{0}}$ and $C^{-}=\frac{\Delta t-2\tau}{2\Delta t\sigma_{0}}$. Replacing $E_{z}^{n}\left(i_{d}+1,j_{d}\right)$ with (23) in (14), and grouping identical components, yields

[TABLE]

with $D^{+}=1+\frac{\Delta t}{\mu_{0}\Delta x}C^{+}$ and $D^{-}=1-\frac{\Delta t}{\mu_{0}\Delta x}C^{-}$. Following the same procedure for $H_{y}^{n+\frac{1}{2}}\left(i_{d},j_{d}\right)$ yields

[TABLE]

As in [14], we set the incident wave as the gaussian-derivative wave $I_{z}(t)=-\sqrt{2e}\beta\left(t-t_{0}\right)\exp\left(-\beta^{2}\left(t-t_{0}^{2}\right)^{2}\right)\>\left(\mu\textrm{A}\right)$ with $\beta=10/t_{0}$, $t_{0}=1\>\textrm{ps}$ and $e$ the electron charge. The position of the source and observation point is shown in Fig. 6.

Simulation results are plotted in Fig. 7 and satisfactorily compared with the result of [14].

VI-D Example 4

The next example is the 2D analysis of the 1D space-time reflection-less metasurface described in Fig. 8. This metasurface is synthesized to linearly oscillate in space and time from fully absorbing to half transmitting as shown in the figure. From the metasurface synthesis equations (1) and Footnote 1, the susceptibilities for an absorbing and half transmitting metasurfaces are $\tilde{\chi}_{mm}^{yy}=\tilde{\chi}_{ee}^{zz}=\frac{-2c_{0}}{j\omega}$ and $\tilde{\chi}_{mm}^{yy}=\tilde{\chi}_{ee}^{zz}=\frac{-2c_{0}}{3j\omega}$, respectively, or generally $-\frac{2c_{0}}{j\omega\alpha(t)}$, where we set the periodic time variation with triangular period

[TABLE]

with $\Delta t=15.349\>\textrm{ps}$. According to this relation, $\alpha(0)=1$, $\alpha(100\Delta t)=3$, $\alpha(200\Delta t)=1$, and so on, periodically. Therefore, the corresponding device is time-varying metasurface linearly oscillating between full power absorption ($\alpha=1$) and half power transmission ($\alpha=3$). The time-varying GSTC for the metasurface is thus

[TABLE]

Following the procedure outlined in Sec. III-B and Sec. IV, the update equation for $H_{y}^{n+\frac{1}{2}}\left(i_{d},j_{d}\right)$ is found as

[TABLE]

with $F^{+}=1+\frac{\Delta tc_{0}}{2\Delta x\alpha(t)}$ and $F^{-}=1-\frac{\Delta tc_{0}}{2\Delta x\alpha(t)}$, while that for $E_{y}^{n}\left(i_{d}+1,j_{d}\right)$ reads

[TABLE]

with $R^{\pm}=1\pm\frac{c_{0}\Delta t}{2\Delta x\alpha(t)}$.

The simulation results are shown in Fig. 9. The incident wave is a normally incident plane wave with gaussian amplitude profile, specifically $E_{z}(t,y)=\sin\left(2\pi ft\right)\exp(-y^{2}/0.02)$ with $f=2\>\textrm{GHz}$. The simulation has been run for 30,000 time steps and no sign of instability has been observed. The transmitted and incident fields right after and before the metasurface, respectively, are shown in Fig. 9. The total field is symmetric with respect to the $y=0$ axis since both the metasurface and the incident field in symmetric. The transmitted field amplitude progressively decreases towards the edges of the metasurface, as expected from its specified spatial profile [Fig. 8]. The total fields in the transmitted ($x>0$) and incident ($x<0$) regions are plotted in Fig. 9. The time variation of the transmitted field may visualized by inspecting the propagation in the plane of maximal intensity, $y=0$; it fully corresponds to the specifications [Fig. 8].

Figure 9 plots the temporal [Fig. 10] and spectral [Fig. 10] profiles of the incident and transmitted fields just before and after the metasurface, respectively. Figure 10 shows that the transmitted field follows the metasurface periodicity ($200\Delta t$). Figure 10 shows the harmonics generated in the transmitted field as a result of time variation, as expected from Sec. II-B.

VI-E Example 5

Finally, we address the 1D problem of the reflection-less half-absorbing and half-transmitting dispersive 0D metasurface. From the previous example, the susceptibilities are $\tilde{\chi}_{\textrm{mm}}^{yy}=\tilde{\chi}_{\textrm{ee}}^{zz}=j\frac{2c_{0}}{3\omega}$. Substituting these susceptibilities into (1) yields

[TABLE]

Following again the procedure outlined in Sec. III-B and Sec. IV, the update equation for $H_{y}^{n+\frac{1}{2}}\left(i_{d}\right)$ is found as

[TABLE]

while that for $E_{z}^{n+1}\left(i_{d}+1\right)$ is

[TABLE]

Results are shown in Fig. 11. As specified, half of the wave is transmitted and the other half is absorbed by the metasurface. Moreover, the zero reflection expected from the equality of the electric and magnetic susceptibilities is also verified.

VII Conclusion

We have proposed an FDTD scheme based on Generalized Sheet Transition Conditions (GSTCs) for the simulation of polychromatic, nonlinear and space-time varying metasurfaces. It has been validated and illustrated by five examples. This scheme is fully numerical and very easy to implement, and represents a fundamental extension of the standard FDTD algorithm to general sheet discontinuities. It is currently being extended by the authors to accommodate general metasurface dispersion.