FDTD Simulation of Dispersive Metasurfaces With Lorentzian Surface Susceptibilities

A Finite-Difference Time-Domain (FDTD) simulation of broadband electromagnetic metasurfaces based on direct incorporation of Generalized Sheet Transition Conditions (GSTCs) into a conventional Yee-cell region has been proposed for arbitrary wave excitations. This is achieved by inserting a zero thickness metasurface inside bulk nodes of the Yee-cell region, giving rise to three distinct cell configurations - Symmetric Cell (SC), Asymmetric Cell (AC) and Tight Asymmetric Cell (TAC). In addition, the metasurface is modelled using electric and magnetic surface susceptibilities exhibiting a broadband Lorentzian response. As a result, the proposed model guarantees a physical and causal response from the metasurface. Several full-wave results are shown and compared with analytical Fourier propagation methods showing excellent results for both 1D and 2D field simulations. It is found that the TAC provides the fastest convergence among the three methods with minimum error.


I. INTRODUCTION
Electromagnetic (EM) metasurfaces are two-dimensional equivalents of volumetric metamaterials and are composed of 2D arrays of sub-wavelength scatterers. By engineering these scatterers across the surface, various impressive waveshaping transformations can be achieved for multiple applications such as generalized refraction, holography, polarization control, imaging, and cloaking, to name a few [1], [2]. Metasurfaces perform such wave transformations as a result of the complex interplay between the electric and magnetic dipolar moments generated by the scatterers, which is sometimes also referred to as a Huygens' configuration [3], [4]. A convenient implementation of such metasurfaces is using all-dielectric resonators, which naturally produce the electric and magnetic dipoles moments, and when properly designed, provide zero back-scattering, resulting in a perfect transmission [5]- [7].
A recently growing area of interest is re-configurable and time-varying metasurfaces, where the surface polarizabilities are real-time tunable. A more general description of such dynamic metasurfaces is a space-time modulated metasurface, where the surface polarizabilities are a function of both The associate editor coordinating the review of this manuscript and approving it for publication was Mohammad Tariqul Islam . space and time, resulting in a traveling-wave type perturbation on the metasurface. They are the 2D equivalents of general space-time modulated media [8], [9], which have found important applications in parametric amplifiers and acousto-optic spectrum analyzers, [10]- [12], for instance. Space-time modulation has led to various exotic effects such as harmonic generation and non-reciprocity [13], [14], that has also been recently explored using metasurfaces [15], [16] for advanced wave-shaping applications. Their attractive features lie in achieving non-reciprocity using purely nonmagnetic materials, which has important practical benefits in engineering systems, related to the high-frequency operation and no requirement of a magnetic bias. Another equally important class is that of nonlinear metasurfaces exhibiting field-dependent susceptibilities for applications in beamshaping and advanced temporal frequency control [17]- [20].
The EM modeling of metasurfaces with such advanced wave manipulations necessitates a need for efficient timedomain simulation of these structures. While practical metasurfaces are sub-wavelength in thickness (δ λ 0 ), they can be efficiently modeled as space-discontinuities, described using electric and magnetic surface susceptibilities, i.e.χ ee (ω) andχ mm (ω), respectively. They thus model practical metasurfaces as zero thickness structures, thereby transforming them into a single-interface problem. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Such space-discontinuities can be rigorously modeled using Generalized Sheet Transition Conditions (GSTCs). Various numerical approaches have been recently presented, where the GSTC conditions are incorporated into bulk Maxwell's equations [21]- [23], including finite-difference formulations in the frequency domain to accurately analyze the transmitted and reflected fields of a general zero-thickness metasurface [24]- [27].
In contrast to frequency-domain modeling, the rigorous time-domain modeling of EM metasurfaces explicitly requires causality considerations for accurate modeling of practical metasurface responses. Since metasurfaces are constructed using sub-wavelength resonators, they are operated around the resonant frequencies where the EM waves have maximum interaction with the metasurface. Consequently, these metasurfaces are naturally very dispersive, i.e. χ e,m (ω) = const. The geometrical shapes of the constituting scatterers are primarily responsible for their resonant behavior, in spite of their design being generally based on non-dispersive materials (typically metals and dielectrics). This operation of the metasurface in a dispersive (and thus broadband) regime, demands a physical description of these resonators consistent with the causality requirements. This in turn, requires a causal description of the equivalent surface χ e,m of the metasurface, in frequency (or time domain), i.e.χ e,m =χ e,m (ω) or χ e,m = χ e,m (t). This requirement is also critical in the accurate time-domain modeling of general space-time modulated metasurfaces and nonlinear surfaces, where new spectral frequency components are generated. This subsequently further necessitates a complete description of the surface χ e,m encompassing these frequencies as well, in addition to the bandwidth of the input excitation.
Consequently, little work has been done in time-domain modeling of metasurfaces, which necessitates these additional requirements of specifying the surface in terms of its temporal dispersion, i.e., frequency-dependent constitutive parameters. The work in [28], [29] presented a GSTC-FDTD formulation of a dispersive metasurface for the first time, where the surface polarizabilities were described using physically motivated Lorentz-Drude models, typical of the physical sub-wavelength unit cell resonators used in constructing practical metasurfaces. Such a specification requires the introduction of auxiliary equations, with the Lorentz-Drude dispersion expressed in the time-domain using second-order differential equations. This approach, known as the Auxiliary Differential Equation (ADE) method, is well understood to be compatible with the standard FDTD update equations [30]. However, the metasurface modeling of [28], [29] did not rigorously integrate the GSTCs inside the Yee-cell, and computes the fields scattering in the transmission and reflection regions separately.
Other FDTD techniques for modeling dispersive metasurfaces have been proposed in works such as [31]- [36]. An approach in some of these works was to restrict the frequency dependence of the surface susceptibilities to simple cases such as χ ∝ 1/ω [32], resulting in limited formulations of the FDTD equations. A more general alternative suggested is the piece-wise linear recursive convolution (PLRC) method [33], [37], which has also been recently applied to model metasurfaces [31]. PLRC is an integral approach, and numerical convolution is needed for its implementation, unlike the ADE method which is based on partial differential equations. However, while methods based on recursive convolution, including PLRCs, provide good accuracy and computational efficiency, the formulation is non-physical, mathematically opaque and requires the use of a partial fraction extraction from a specified frequency response -which can pose problems with respect to passivity and causality. The method is also limited to linear surfaces.
Conversely, the ADE method has a number of natural advantages for modeling metasurfaces. For a complex surface response, the extraction of the multiple Lorentz-Drude resonances is relatively simple due to the underlying physics being captured by the approach [21]. As the individual Lorentz-Drude response is causal and, if required, passive, the total response is guaranteed to be physical. The method is also attractive as it can naturally be extended to model nonlinear effects and thus is also a more suitable choice to simulate space-time modulated metamaterials and metasurfaces [38]. Moreover, the ADE method has identical accuracy and memory requirements for Lorentz media as the PLRC [39], [40].
In this work, we develop a rigorous Finite-Difference Time-Domain (FDTD) method, where the GSTCs are directly integrated into the FDTD Yee-cell. Compared to the explicit FDTD model presented in [28], [29], [38], where the metasurface is treated as a boundary of a given simulation domain, the proposed method treats the metasurface as an EM scattering entity and is able to process arbitrary broadband excitations. In this context, various strategies for integrating GSTCs in a conventional FDTD Yee-cell are proposed and compared here, assuming Lorentzian surface polarizabilities, which are naturally causal and rigorously captures the fundamentally dispersive nature of typical EM metasurfaces. Moreover, issues related to the exact placement of the metasurface inside the Yee-cell relative to the standard field node, are comprehensively discussed, which have important implications in achieving desired computational accuracy.
The paper is structured as follows. Section II presents the GSTC formulation of zero thickness metasurfaces, and establishes analytical models for specific cases, for benchmarking purposes. Section III proposes three possible different Yee-cell configurations where GSTCs can be incorporated in conventional cells. Section IV shows several simulated results corresponding to these configurations, providing detailed comparisons between them. Finally, Sec. V provides concluding discussions and remarks on the applicability of the Lorentzian surface susceptibilities, and conclusions are provided in Sec. VI. Some important field derivations of the proposed Yee-cell configurations are also provided in the Appendix.

FIGURE 1.
A typical configuration of a zero-thickness uniform Huygens' metasurface located at z = 0. It consists of orthogonal electric (p) and magnetic (m) dipole moments, excited with a normally incident plane-wave, resulting in reflected and transmitted fields governed by (1). For simplicity, normal polarization is assumed without any variation of the fields along the y −axis.

II. METASURFACE DESCRIPTION A. GENERALIZED SHEET TRANSITION CONDITIONS (GSTCs)
A zero thickness metasurface, such as the one in Fig. 1, is a space-discontinuity. The rigorous modeling of such discontinuities based on Generalized Sheet Transition Conditions (GSTCs) was developed by Idemen in [41], which were later applied to metasurface problems in [42]. The modified Maxwell-Faraday and Maxwell-Ampere equations can be written in the time-domain as, where ψ represents the differences between the fields on the two sides of the metasurface for all the vector components of the field ψ, i.e., H or E fields. The other terms P s and M s represent the electric and magnetic surface polarization densities, in the plane of the metasurface, which depend on the total average fields around the metasurface [43]. They are defined by,P whereχ ee andχ mm are the frequency dependent electric and magnetic susceptibilities andẼ s andH s are the average EM fields at the surface. The EM coupling related to the bi-anisotropic term is assumed zero here, for simplicity. Furthermore, the surface susceptibilitiesχ ee andχ mm are also treated a scalars for simplicity, as opposed to their most general tensorial forms that account for more general wave transformations.

B. SURFACE POLARIZATION DENSITIES
A primary concern in modeling the metasurface response is that a physical representation of the surface polarizations must be consistent with causality. The metasurface unit cells will have a number of natural resonances, and this response must be captured for the correct broadband response to be predicted. These resonances are naturally modeled by Lorentzian functions, and a summation of correctly parameterized Lorentzians is an appropriate model of the surface, chosen in this work, i.e.
A key consideration in the use of Lorentzians is that they represent a physical process and therefore are implicitly causal. Besides, they naturally take into account the dispersive effects of the metasurface, which have practical importance in the EM interaction of metasurfaces with broadband excitations. While a typical metasurface requires several Lorentzian contributions to accurately model broadband surface susceptibilities, such as for all-dielectric unit cells [28], here we assume a single resonant Lorentz response for the electric and magnetic polarizations, for the sake of clarity and simpler forthcoming analytical expressions, so that where ω p , ω 0 , and α are the plasma frequency, resonant frequency, and the loss-factor of the oscillator, respectively. The subscripts e and m denote electric and magnetic quantities. It should be noted that the use of a single Lorentzian here is strictly for simplicity. Therefore, it is possible to extend the proposed method for multiple resonance contributions due to linear field superposition.
To formulate a time-domain Yee surface cell, we need a time-domain representation for the surface polarizations and must transform (4) to the time domain, using the inverse Fourier transform, leading to It is convenient to formulate these two equations as a 1 st order system using two variables. For example, for the case of Fig. 1, defining P s and ω e0 P s = dP s /dt along the y−axis, we obtain, allowing us to write the first equation as, In a similar manner we can define M s and ω m0 M s = dM s /dt and obtain, Let us first consider a specific case of a linear and uniform metasurface, which is excited with a normally incident pulsed plane-wave. It represents a simple case whose time-domain transmitted, and reflected fields can be obtained using a standard Fourier propagation method [44]; which provides a baseline for comparison to the Yee cell simulations. Now, consider a Huygens' metasurface illuminated with a normally incident plane-wave, as shown in Fig. 1. For simplicity, but without loss of generality, the problem is assumed to be 2D, where all y−variations are considered to be zero. Assume a monochromatic excitation with a frequency ω and expressingẼ s andH s in terms ofẼ t (ω),Ẽ r (ω) and E 0 (ω); the frequency domain representation of the transmitted, reflected and incident plane-waves in the metasurface, respectively. Substituting these corresponding fields into (2) leads to,P For a normally incident plane-wave, (1) becomes Equations (9a) and (9b) can further be placed into a matrix form by defining, Using this equation after specifyingẼ 0 (ω) = F{E 0 (0 − , t)}, we can then determine E t (0 + , t) and E r (0 − , t) using an inverse Fourier transform, i.e.
which finally represents the time-domain reflection and transmission of the metasurface.

III. YEE CELL FORMULATION
In this section, we will describe the integration of GSTCs into bulk Yee-cells, using three possible configurations. The formulation will be developed for 2D propagation in the x − z plane. However, the proposed approach can be straightforwardly extended into a full 3D simulation. The entire simulation region consists of two types of nodes: bulk nodes for modeling the reflection and transmission region following Maxwell's equations, and surface-nodes modeling the zero thickness metasurface following the GSTCs. For a bulk 2D Yee-cell defined for propagation on the x − z plane (see Fig. 2), the basic equations for the electric and magnetic fields, are given by conventional update equations, where the subscript is the spatial position (with i and j representing node positions in z and x respectively), and the superscript denotes the time-step. A field evolution is then obtained by stepping in time using the following procedure: 1) Update the H's at n + 1/2 (t n+1/2 = t n + t/2) using previous E's and the boundary conditions. 2) Update E's at n + 1 (t n+1 = t n + t) using the H's calculated at n + 1/2. A key feature of this procedure is that it determines the unknown current fields using only previously calculated fields. For example E n+1 's are determined from the H n+1/2 's. Therefore the method is a strict explicit time marching method. As an explicit method the spatial step, , and the time step, t, should always satisfy the Courant condition u where u is the local speed of light) to guarantee stability [30].
Next is the metasurface region. When forming the update equations for the metasurface cells, the nature of the GSTC equations in (1) keep the H x | n+1/2 and the E y | n+1 fields coupled, which not allows them to be solved in a simple sequential manner. This is due to the polarizations being naturally solved at the time point associated with their fields. Thus, we solve for M s | n+1/2 and P s | n+1 by formulating a selfconsistent solution to the unknowns present in each cell and then update all fields after a complete time step.

A. YEE-CELL CONFIGURATIONS
In this section, we present three different possibilities to incorporate a zero thickness surface via GSTCs inside a bulk Yee-cell region with minimal disruption. The three following surface cells we wish to consider are shown in Fig. 2: 1) Symmetrical Surface Cell (SC): This cell is identical in form to the cell used in [28] for formulating an explicit surface as an internal boundary condition. The surface is inserted midway between two electric field nodes (E y | k−1 and E y | k ) and the H x node present at the position k-1/2 is split into two nodes on either side of the surface denoted as H x | s − and H x | s + . 2) Asymmetrical Surface Cell (AC): The second cell is formed by simply inserting the surface midway between an H x node and E y node. This requires no new nodes to be defined but produces an asymmetrical cell structure.

3) Tight Asymmetrical Surface Cell (TAC):
The third structure is similar to the AC cell but inserts new nodes either side of the surface. On the left side a new electric field node is inserted (E y | s − ) and on the right a magnetic field node is inserted (H x | s + ).

B. YEE-CELL UNKNOWNS
A self-consistent solution for the surface cell requires a clear identification of the variables defined at and near the surface.
Although most of the nodes are essentially bulk nodes, there are still irregular nodes that require a special update equation. See, for example, the E y | k node in the SC cell, it relies on the value of H x | s + ,j that is not known until the surface cell is solved. Assuming all bulk nodes (H x | n+1/2 , H z | n+1/2 and E y | n+1 ) were updated, we obtain the following unknown vectors, X, by inspecting each of the surface cells: Given these unknowns, we now need to formulate a linear set of equations with the GSTCs of (1) and the polarization responses in (6) and (7), so they can be integrated into the Yee cell algorithm. The update equations for the surface cells derived in the following sections are complex requiring the self-consistent solution of variables at two times (n + 1/2 and n + 1). To clarify them we have used boxed variables such as E y | n+1 k,j to denote the unknowns. Other values will be known at the time of solution for the fields in the cell.

C. METASURFACE FIELD EQUATIONS
The first GSTC equation (1a) specifies a relationship between the field across the surface ( E y ) and the magnetic polarization. To discretize this equation we have two choices: (a) As a first choice, we naturally impose the electric field across the surface at the time point n + 1 and use central difference in time for the polarization M s (which will be solved on the half steps n + 1/2) providing, This equation is problematic, however, as M s | n+3/2 j is in the future and not part of our solution. We must therefore make an approximation and use, (b) The second choice is to center E y at the time step n, allowing us to write, However, this choice is also not useful, as for small χ ee , the above equation becomes problematic. This can be seen by setting χ ee = 0, which produces, E y | n s = 0. Such an equation does not provide a relationship between any of the unknowns present in E y | s and will produce an VOLUME 8, 2020 under-determined set of cell equations. We will therefore use (11) in the subsequent formulation. The use of this assumption may require the spatial and time steps of the simulation to be smaller then that typically used for Yee based FDTD modeling and as we will see in the results section where convergence with respect to discretization is investigated the consequence of modeling the dispersive surface is a need to use relatively fine meshes.
The first GSTC equation (1a) relating H x to the electric polarization can be handled more straightforwardly. Using the H x values at the surface centered on the time step n + 1/2, we have, Both of these equations (11) and (12) where the difference is obtained by simply using the appropriate nodes on either side of the surface.

D. SURFACE POLARIZATION EQUATIONS
The time domain surface polarization equations (6) and (7) need to be discretized in time. As these equations represent a 2 nd order system, a trapezoidal formulation was used producing, For the three cells these two equations only differ in the nature of the forcing terms E y | s,j and H x | s,j -the average field at the surface. For the three cells, we have by inspection, It is important to note that for the two asymmetrical cells the fields applied to the surface E y | s,j and H x | s,j are not coincident in space. For the AC cell the E y | s,j is applied at the position k-1/2 and H x | s,j at k. For the TAC cell E y | s,j is applied at the position k-1/4 and H x | s,j at k-3/4. This mismatch can be expected to produce some error in the It is largest for the AC cell and a possible modification is to use, which brings the forcing function into alignment at k − 1/2. Unless otherwise noted, the AC cell will use the symmetrical forcing formulation.

E. SPECIAL CELL UPDATE EQUATIONS
In addition to equations defining the GSTCs and the polarizations, each cell has a number of nearby nodes which, although placed in the bulk, are dependent on the surface equations and thus must be solved at the same time.
(a) For the SC, we have special update equations for the nodes E y | n+1 k−1,j and E y | n+1 k,j , given by The AC cell only requires one special update for E y | n+1 k−1,j , given by For each cell configuration, the equations (11-13) and one of (17), (23) or (19) form a complete set of linear equations defining the surface variables. These equations can be assembled into a compact matrix form, given by For the specific case of the TAC cell, [ s ] and [F s ] are presented in the appendix as an example. A field evolution is then simply obtained by stepping in time using the following procedure: (a) Update the bulk H's at n+1/2 (t n+1/2 = t n + t/2) using previous E's and the boundary conditions. (b) Update bulk E's at n + 1 (t n+1 = t n + t) using the H's calculated at n + 1/2. (c) Update surface H s 's, and M s at n + 1/2 (t n+1/2 = t n + t/2) and E s 's and P s s at n + 1 (t n+1 = t n + t) using

IV. VALIDATION
To test the various Yee-cell configurations, the algorithm above was integrated into a standard 1D/2D Yee-cell based simulator, with a metasurface configured as an SC, AC or TAC cell. The simulation setup was configured as a source applied to the left side of the computation region, and perfectly matched layers (PMLs) on the other three sides. A metasurface was placed vertically half way along the z-axis, at z = 0. The intent of this section is to evaluate the robustness and accuracy of the three cell configurations, and do a detailed comparison between them.

A. TRANSPARENT SURFACE
To provide an initial evaluation of the surface cells, simulations were performed for transparent surfaces (or no metasurface), where P s = M s = 0. For such a case the GSTCs reduce to which gives, As of course, the polarization equations are no longer relevant for a transparent region with no metasurface, which can now be represented by one of the update equations from (17)(18)(19) and equation (21). More specifically, they are given for each of the three Yee-cell configurations as follows: (a) For the SC cell, we obtain, These are identical to bulk equations in absence of a metasurface, however, the relationship E y | n+1 k,j = E y | n+1 k−1,j imposed by the GSTCs (21) is not correct for this configuration, and we can thus expect some errors to be present. (b) For the AC cell we get, This is again the bulk update equation. Unlike the SC case, the relationship E y | n+1 k,j = E y | n+1 k−1,j produces the correct update equation this time for the adjacent nodes. The surface is essentially removed from the simulation VOLUME 8, 2020 mesh producing the bulk update equations, and we can expect perfect transparency. (c) For the TAC cell, we obtain the following update equations, These are identical to the bulk update equations except for a 0.75 factor due to a shortening of the cells either side of the surface. This slight distortion of the cell is due to the use of complimentary E and H field nodes on either side of the surface. We can expect near perfect transparency from this configuration. Figure 3(a) shows the computed fields in the reflection (z < 0) and transmission region (z > 0), when the metasurface is numerically removed (χ ee (ω) =χ mm (ω) = 0). The figure shows the response of all three Yee-cell configurations (SC, AC and TAC), and compared with analytical Fourier propagation method of Sec. II-C. The modulation frequency of the input pulse is 230 THz throughout this paper. The transmitted pulse shows a similar response for all Yeecell simulations with a slight discrepancy from the Fourier propagation result. The spatial step size was = λ 0 /100 and was chosen to produce a negligible amount of numerical dispersion. As this dispersion is the same for all cases, this slight pulse distortion can be attributed to bulk effects, which can be reduced by lowering the spatial step size. The reflection, of course, should ideally be zero, and is found to be negligible for AC and TAC cases. However, the SC exhibits a significant amount of reflection due to not producing the bulk update equations correctly (see (22)).
Next, a stepped continuous-wave (CW) Gaussian beam with a width of λ 0 was launched from the left side, with strongly divergent wavefronts. Figure 3(b) shows the computed steady-state E-field distribution obtained for each of the three Yee-cell configurations. As can be seen, no significant reflection is present except in the case of a SC configuration. We can conclude from these observations, that for transparent or nearly transparent surfaces (i.e, χ ee, mm ≈ 0), the SC is inappropriate due to existence of spurious numerical reflections.

B. LORENTZIAN METASURFACE
Let us now introduce a metasurface inside the computational region, at z = 0. Fig. 4(a) shows the transmitted and reflected pulses, corresponding to a normally incident Gaussian pulsed plane-wave on a matched surface (χ ee =χ mm ) for three As can be seen all three surfaces produce similar results. The pulse is strongly dispersed by the presence of a metasurface with Lorentzian susceptibilities. Furthermore, the metasurface response matches very well with the FT propagation result albeit with some numerical dispersion determined by the spatial step size. A matched metasurface should produce no reflections, and the FT result does show this. For the larger step sizes all of the surfaces produce non-negligible reflection with the AC cell, arguably producing the most and the TAC  Fig. 3. The metasurface susceptibilities are defined using a single Lorentzian, for simplicity, with the following parameters: ω ep = 3.01 × 10 11 rad/s, ω e0 = 2π (230 THz) and α e = 7.54 × 10 12 . For the mismatched metasurfaces ω e0 − ω m0 = 2π (15 GHz). The input pulse is Gaussian with full-width-half-maximum (FWHM) of 10 −15 s with a modulation frequency of 230 THz.
the least. For all the Yee-cell configurations, the field reflection can be reduced to a negligible value by decreasing the spatial step size to = λ 0 /400. Figure 4(b) further shows the transmitted and reflected fields for a mismatched metasurface (χ ee =χ mm ) with a spatial step size of = λ 0 /400, where a significant amount of reflection is expected. All of the Yee-cell configurations produce an excellent match to the FT result, with again the TAC cell producing the least amount of error. Finally, to further evaluate the impact of the forcing function in the Lorentzian susceptibilities, a simulation using the AC cell with asymmetrical forcing (15) is also shown in Fig. 5. It can be clearly seen that even for a very small step size ( = λ 0 /400), significant reflections are produced. We can therefore conclude that Lorentzian polarizabilities must be excited with symmetrical forcing functions for all three Yee-cell configurations. Figure 6 further shows the impact of the spatial step size on the scattered fields for both cases of matched and mismatched metasurfaces. The total normalized EM energy is computed in both transmission and reflection regions as a function of step size , and compared with the ideal results obtained using the Fourier transform propagation method. A monotonic convergence is observed for all the three Yeecell configurations, and in particular, AC and TAC are seen to converge faster than the SC, in all cases. It should be noted that the total normalized energy in the computation region is less than 1, due to non-zero α's accounting for dissipation losses in the metasurface.  Fig. 4, computed using the asymmetric cell (AC), using a symmetric and asymmetric forcing function in the surface polarization.
Finally, Fig. 7 shows the scattered fields from a matched and mismatched metasurface, excited with a stepped CW Gaussian beam, computed using the TAC configuration. As expected, the matched surface produces a phase discontinuity at the surface but with no disruption of the beam propagation. On the other hand, the mismatched metasurface creates a more complex response as the reflected waves from the surface cause a standing wave to form in the reflection VOLUME 8, 2020  region, as compared to purely forward propagating waves in the transmission region.
While the use of the Lorentzian metasurface guarantees a causal physical system, the numerical technique employed in the description of the surface could introduce numerical instabilities. Detailed simulations of the case of Fig. 7 were run to harmonic steady state for up to 0.25 ns (roughly 6000 optical cycles) without observing any instabilities. Moreover, for broadband temporal pulses, such as the ones investigated here, this approach has been found to be adequate.

V. CAUSALITY CONSIDERATIONS FOR CONSTANT SURFACE SUSCEPTIBILITIES
The prime motivation for modeling metasurface susceptibility densities using Lorentzian profiles inside of a Yee-cell region, is to incorporate a physically motivated causal response from the metasurface sub-wavelength unit cells. While an assumption of a constant susceptibility will greatly simplify the proposed Yee-cell algorithm, it fails to rigorously capture two important physical effects: a) causality, and b) dispersion. While it is clear that constant surface susceptibilities are inherently non-dispersive, the causality aspect is not straightforward. This aspect is clarified in this section, to emphasize the importance of using Lorentzian susceptibilities in the proposed FDTD method.
Let us consider a lossless uniform matched metasurface [χ ee (ω) =χ mm (ω) = χ 0 , where χ 0 ∈ R] excited with a pulsed Gaussian plane-wave, whose envelope is given by The transmission function of a matched metasurface is given by (9), as where the last equality is based on the fact that |T (ω)| = 1 ∀ ω, i.e., an all-pass transfer function. For small arguments of the inverse tangent function, T (ω) ≈ exp{−jωχ 0 /c}, so that the output of the metasurface is given by where k 0 = ω/c. Therefore, the metasurface output is also a Gaussian pulse, as expected. However, its peak is now located at a time instant t peak = k 0 χ 0 . There are now two following possibilities: 1) χ 0 > 0: The output pulse is located at t peak ≥ 0, i.e., a positive time-delay. 2) χ 0 < 0: The output pulse is now located at t peak < 0, i.e., a negative time-delay or a time advance. While Case 1 is naturally causal, Case 2 represents a non-causal response, as the output pulse appears before the input. Therefore, for this simple case of a matched metasurface, it can be concluded that negative and constant surface susceptibilities, represent a non-physical system, and thus are not allowed. This is consistent with the fact that causal EM metamaterials with negative constitutive parameters < 0 and µ < 0, must be dispersive to allow positive time-average stored electric and magnetic energies [45], [46].
To numerically demonstrate this, consider a broadband pulsed excitation of a uniform metasurface with a very short Gaussian pulse. Fig. 8 shows the computed response corresponding to both analytical Fourier transform propagation approach and the proposed Yee-cell using the TAC configuration, for the two cases when Re{χ 0 } > 0 and Re{χ 0 } < 0. As expected, the TAC simulation and the FT simulation match perfectly for the first case providing a positive time delay. The second case with negative real part of χ 0 produces a phase advance for the transmitted field, thereby creating a non-causal response in the FT simulation with the transmitted pulse running ahead of the excitation envelope. As we would expect, for the Yee-cell simulation, this causes an instability which is indicated by the amplification and gross distortion of the transmitted pulse [47]. This further demonstrates the need for the use of a physical causal surface representation, such as a Lorentzian response in a numerical Yee-cell model.

VI. COMPUTATIONAL IMPACT
The presence of a metasurface in the simulations above has an impact on the simulation time due to an increase in computational intensity. The surface Yee-cell updates are considerably more numerically intensive then a simple bulk Yee-cell update. However, this is mitigated by the fact that there are far fewer of the surface cells than the bulk cells. A primary consideration is that for simple geometries the number of surface cells will scale linearly with the resolution of the meshing, whereas the number of bulk cells will increase geometrically.
This trade-off is illustrated in Fig. 9 where for two geometries we present the simulation time as a function of the spatial step size. The first geometry is a quasi-1D simulation of 20 µm long strip with a fixed number of cells in the y direction, and variable number in the x direction. The second geometry is a 2D simulation of a rectangular region 20 × 20 µm 2 as shown in Fig. 4. For all simulations the physical dimensions were kept constant, and simulations were run for increasing mesh densities from 12 to 95 steps per wavelength for all surfaces and the case of no surface, using the later as a reference case. As can be seen from the figure, for the quasi-1D case the presence of a surface produces a noticeable increase in simulation time for low mesh densities. In this case the computation overhead is constant as the number of surface cells is fixed.
For the second geometry, the situation is similar, but more pronounced. At low mesh densities the CPU cost of the surface calculations is substantial, however, as the mesh density increases, the number of bulk cells increases geometrically, as opposed to only a linear increase of number of surface cells. Thus, for moderate or fine meshing, the CPU cost of the surface becomes almost insignificant. In both examples, the computation overload of the surface, expressed as a percentage of the total simulation time, decreases as the mesh density increases, due to larger number of bulk cells.

VII. CONCLUSIONS
An FDTD simulation of broadband electromagnetic metasurfaces has been proposed based on direct incorporation of GSTCs inside of Yee-cell region, for arbitrary wave excitations. This has been achieved by inserting a zero thickness metasurface inside bulk nodes of the Yee-cell region, giving rise to three distinct cell configurations -SC, AC and TAC. In addition, the metasurface has been modelled using electric and magnetic surface susceptibilities exhibiting a broadband Lorentzian response. As a result, the proposed model guarantees a physical and causal response from the metasurface. Several full-wave results have been presented, and compared with analytical Fourier propagation methods showing excellent results, for both 1D and 2D fields simulations. It has been further found that the TAC provides the fastest convergence among the three methods with minimum error. While the case of scalar surface susceptibilities exhibiting a single Lorentzian resonant contribution is described here, and was found to be adequate, the proposed method can be extended to full tensorial susceptibilities with an arbitrary number of Lorentzians within a fully 3D simulation environment. Finally, the Lorentzian description surface polarizabilities enables a physical, intuitive and practically useful description of the metasurface, thereby easily extendable to space-time modulated metasurface cases, where for instance the resonant frequency and the loss-coefficient of the Lorentz function can be easily modulated in time [38]. In such cases, the ADE-FDTD technique becomes a preferred choice compared to other FDTD techniques for modeling material dispersion, including PLRC methods. Future works include extension of the proposed methods to model more exotic metasurfaces and performing its detailed stability analysis.

A. TIGHT ASYMMETRICAL SURFACE CELL MATRIX FORMULATION
To illustrate the formation of field matrix equation (20), let us take an example of a tight asymmetrical cell among the possible three configurations of Fig. 2. The equations describing the tight asymmetric cell are given by (11)(12)(13) and (23). Using them, we obtain the following set of field equations: 1) The pair of GSTCs for both E-and H -fields.   These equations can be placed in appropriate matrix forms as given in (24)- (26), as shown at the bottom of the previous page, which can be now used to update the subsequent fields on the Yee-cell nodes.