Efficient Dispersive GSTC-FDTD Algorithm Using the Drude Dispersion Model

Metasurfaces are artificial sheets with sub-wavelength thickness and they are two-dimensional equivalents of metamaterials. The generalized sheet transition conditions (GSTCs) have been recently proposed for electromagnetic analysis of the metasurfaces. In GSTCs, the metasurface is generally modeled as a sheet with zero-thickness. However, the conventional finite-difference time-domain (FDTD) method is not straightforwardly applied to analyze electromagnetic wave propagation in the metsurface by harnessing GSTCs because GSTCs exhibit electric and magnetic discontinuities. Alternatively, the GSTC-FDTD formulation is highly suitable for analyzing the electromagnetic properties of metasurfaces by introducing electric and magnetic virtual grids. Meanwhile, metasurfaces can be realized by using 2-D materials such as black phosphorus and thus the dispersion characteristics of metasurfaces should be considered. In this work, we propose an efficient dispersive GSTC-FDTD algorithm by employing the Drude dispersion model. Moreover, for the first time, the numerical surface susceptibility inherent to the dispersive GSTC-FDTD formulation is derived and its numerical accuracy is investigated. Numerical examples illustrate high efficiency of the proposed Drude-dispersive GSTC-FDTD algorithm.


I. INTRODUCTION
There are various computational methods for analyzing complex electromagnetic problems [1]- [8]. Among many numerical methods, the finite-difference time-domain (FDTD) method has been popularly employed to analyze electromagnetic wave propagation in complex media such as lossy, anisotropic, non-linear, or dispersive media [9]- [18]. The FDTD method is a grid-based numerical analysis technique that discretizes spatial and temporal partial derivatives of the time-dependent Maxwell's curl equations using the central-difference scheme (CDS) [19], [20]. The remarkable advantages of the FDTD method are robustness, simplicity, and broadband analysis in a single simulation.
The metasurfaces are sub-wavelength-thickness sheets to control electromagnetic wave propagation to yield extraordinary properties. Metasurfaces have been successfully applied to a broad variety of applications, such as flat-lenses [21], [22], hologram generation [23], remote processing [24], The associate editor coordinating the review of this manuscript and approving it for publication was Debdeep Sarkar . harmonic generation [25], dual-band antennas [26], and polarization analyzers [27]. The thickness of metasurfaces is much thinner than a wavelength in general and thus they can be modeled as zero-thickness sheets. The generalized sheet transition conditions (GSTCs) was proposed to efficiently model these metasurfaces using the surface susceptibility [28]. The conventional FDTD method cannot be simply used to analyze metasurfaces by harnessing the efficient GSTCs since the standard E-and H -FDTD grids are not located on zero-thickness sheets.
Recently, the GSTC-FDTD algorithm was proposed by introducing virtual E-and H -FDTD grids and considering the surface susceptibility between virtual FDTD grids [29], [30]. Vahabzadeh et al. analyzed polychromatic, nonlinear, and space-time varying metasurfaces by using the GSTC-FDTD algorithm [30]. In addition, the GSTC-FDTD algorithm was extended to analyze electromagnetic properties of dispersive metasurfaces by adopting the Lorentz dispersion model [31]. However, this dispersive GSTC-FDTD formulation leaves some regrets in terms of computational efficiency.
In order to improve computational efficiency, we propose a GSTC-FDTD formulation based on the Drude dispersion model. Our proposed GSTC-FDTD formulation can accurately analyze the electromagnetic characteristics of dispersive metasurfaces, with high computational efficiency. The remainder of the paper is organized as follows. In the next section, the GSTC-FDTD algorithm is briefly reviewed and then the GSTC-FDTD formulation by applying the auxiliary differential equation (ADE) to the Drude dispersion model is proposed. Next, we investigate the numerical surface susceptibility of the proposed Drude-based GSTC-FDTD and the previous Lorentz-based GSTC-FDTD formulations. Numerical examples are employed to illustrate superiority of the proposed Drude-dispersive GSTC-FDTD formulation over the previous Lorentz-dispersive GSTC-FDTD formulation and the conventional FDTD algorithm. Finally, concluding remarks are followed.

II. PROPOSED DISPERSIVE GSTC-FDTD FORMULATION
Let us consider a metasurface located on the xy-plane of the Cartesian coordinate at z = , as shown in Fig. 1. In the GSTC-FDTD algorithm, a virtual electric FDTD grid located at − and a virtual magnetic FDTD grid located at + are introduced. The GSTC-FDTD update equations for E x and H y field components are expressed as follows where ε 0 and µ 0 indicate the permittivity and the permeability of the free space respectively. Note that z indicates the FDTD grid size along the z direction, k indicates the space index along the z direction, t denotes the FDTD time step size, and the superscript indicates the time index. Here, the virtual electric field, E x ( − ), and the virtual magnetic field, H y ( + ), are calculated through the GSTCs algorithm. In the frequency domain, the GSTC algorithm can be written as (an e jωt time dependence is assumed) [28] − H y (ω) = J xx ee (ω) + J xy em (ω), where J xy em (ω) = jk 0 χ xy em (ω) H y,av (ω), M yy mm (ω) = jωµ 0 χ yy mm (ω) H y,av (ω), Note that ω is the angular frequency, k 0 is the propagation constant of the free space, χ represents the frequency-dependent surface susceptibility tensor, with ϕ representing the spectral electric or magnetic fields, and the superscripts ''inc'', ''tra'', and ''ref'' denote the incident, transmitted, and reflected waves, respectively [32]. In the FDTD framework, it is necessary to convert (3) and (4) into the time-domain counterparts. By applying the inverse Fourier transform (IFT) and the central averaging scheme, we obtain the following update equations for the GSTCs algorithm in the discrete FDTD world as follows Substituting (9) into (1), and (10) into (2) and then rearranging them, we have the following GSTC-FDTD update equations for E x and H y : Now, let us consider the dispersive properties of the metasurfaces. In this work, we employ the Drude dispersion model to derive an efficient dispersive GSTC-FDTD formulation. We show details on how to derive update equations for the Drude-dispersive GSTC-FDTD algorithm by considering only the electrical polarization current density, J xx ee (ω). Using VOLUME 10, 2022 the Drude dispersion model for the surface electric susceptibility, J xx ee in (5) can be expressed as follows By rearranging (13) and applying the IFT, we have Note that the jω terms in the numerator and the denominator in (13) are cancelled out each other and thus the first-order time derivative is involved in the ADE for the Drude dispersion model. By applying the standard FDTD discretization to the above equation, we have Rearranging (15), we obtain the following ADE update equation for J xx ee : Now, by substituting (16) and (17) into (11) and (12) respectively, we have Here, the average electric and magnetic field components are defined as E n+1  (18) and (19) and rearranging them, we finally obtain the proposed Drudedispersive GSTC-FDTD update equations for E x and H y : On the other hand, as alluded previously, the Lorentzdispersive GSTC-FDTD formulation was introduced in [31]. The surface susceptibility of the Lorentz dispersion model is expressed as It is worthy to include the update equations for the Lorentzdispersive GSTC-FDTD formulation in [31] and they are as follows It should be noted that ω 2 0,ee is involved in the Lorentz dispersion model and thus the second-order time derivative should be considered in Lorentz-dispersion-based polarization current densities, which leads to worse computational efficiency versus the Drude dispersion model counterpart. Table 1 shows the number of field variables and operation count for each dispersive GSTC-FDTD formulation. As shown in Table 1, the computational efficiency of the proposed Drude-dispersive GSTC-FDTD formulation is better than in the Lorentz-based counterpart.

III. NUMERICAL SURFACE SUSCEPTIBILITY
In this section, the numerical surface susceptibility is derived for the proposed Drude-dispersive GSTC-FDTD algorithm and the Lorentz-dispersive GSTC-FDTD algorithm. Since the FDTD method has the discrete nature, the numerical surface susceptibility of the media is different from the analytical surface susceptibility.
First, we consider the numerical surface susceptibility for the proposed Drude-dispersive GSTC-FDTD formulation. For simplicity, we derive the numerical surface susceptibility for the electrical polarization current density only. Let us express E n = Ee jωn t and J n = Je jωn t [12], [33], [34].
By rearranging (28), we obtain In above, ω = 2 t tan ω t 2 . We finally derive the numerical surface susceptibility of the proposed Drude-dispersive GSTC-FDTD algorithm: Next, let us derive the numerical surface susceptibility of the Lorentz-dispersive GSTC-FDTD formulation. In the similar fashion, we have By dividing both sides by e jωn t and rearranging the resulting equation, we obtain where ω 2 0,ee = ω 2 0,ee /(cos ω t 2 ) 2 . Finally, we obtain the numerical surface susceptibility of the Lorentz-dispersive VOLUME 10, 2022
In [40], Feng et al. proposed BP-based metasurfaces to achieve near-unity infrared absorption by controlling the angle of incident light and polarization. In [41], the authors employed BP as metasurfaces exhibiting broadband and polarization-insensitive coherent perfect absorption in THz band. BP is a promising 2-D material for metasurface applications because it can yield extraordinary electromagnetic properties such as tunable anisotropic plasmonic responses [42]. In this work, as a proof of concept, we consider BP to validate the numerical accuracy and efficiency of the proposed Drude-dispersive GSTC-FDTD algorithm. The anisotropic surface conductivity of BP can be calculated by employing a semiclassical Drude model [35] as where j = x, y represents each in-plane direction. In above, η e is the relaxation rate, D j = πe 2 n sj /m j is the Drude weight, n sj is the electron doping, m j is the electron mass along the j direction, and the electron mass along the x direction and y direction can be evaluated by whereh is the reduced Plank constant, γ = 4d/π , η c = h 2 /(0.4m 0 ), υ c =h 2 /(1.4m 0 ), BP is the thickness dependent direct band gap, d is the scale length of BP, and m 0 is the static electron mass. In this work, BP = 2 eV, η e = 10 meV, n sj = 10 13 /cm 2 , and d is 0.223 nm [35]. The electric surface susceptibility of BP can be expressed in the Drude dispersion model as follows where ω D,ee = D j ε 0 π and γ ee = η e /h. First, we investigate the numerical surface susceptibility of the proposed Drude-dispersive GSTC-FDTD algorithm and the Lorentz-dispersive GSTC-FDTD algorithm in the frequency range of 1 THz to 10 THz. Fig. 2 shows the analytical surface susceptibility and its numerical counterparts for the proposed Drude and the conventional Lorentz dispersive GSTC-FDTD algorithms. Note that the FDTD grid size was set to 1 nm in 2-D materials in literature [13], [43] but, in this work, we set the FDTD grid size as 10 nm (which is 3000 cells per minimum wavelength in free space) for the dispersive GSTC-FDTD algorithm, unless specified otherwise.
The time step size is set to 3.3 × 10 −17 s for the CFL stability condition [19]. The numerical susceptibilities are in good agreement with the analytical surface susceptibility. Fig. 3 shows the relative errors of the numerical surface susceptibility for the Drude-and Lorentz dispersive GSTC-FDTDS and it is observed that both errors are same and they can be negligible.
Next, we analyze electromagnetic properties of BP using actual FDTD simulations. In FDTD simulations, the one-dimensional computational domain has a physical length of 100 nm along the z direction and the metasurface sheet (BP) is located at the center of the computational domain. The electric current with a differentiated Gaussian pulse in the time domain is used for source excitation. The whole FDTD domain consists of 30 cells including 10-cell perfectly matched layers (PMLs) [44]- [46] on both ends. All FDTD simulations are performed on Intel i7-10700 CPU. Let us consider the BP sheet when the direction of the electron mass is armchair or zigzag. When the direction of the electron mass is armchair, the x-polarized electric current source is excited and when the direction of the electron mass is zigzag, the y-polarized electric current source is excited. Other simulation setup parameters are maintained. Figs. 4, 5, and 6 shows the reflection coefficient and the transmission coefficient of the BP sheet when the FDTD grid size is 1 nm, 2 nm, and 10 nm. As shown in Fig. 4, the numerical results of the Lorentz-dispersive GSTC-FDTD, the Drude-dispersive GSTC-FDTD, and the conventional Drude-ADE FDTD [19] agree very well with the analytical values [32] for both armchair and zigzag directions when the FDTD grid size is 1 nm. However, when the FDTD grid size is increased to 2 nm or 10 nm, the conventional FDTD results are significantly different from the analytical values, whereas both dispersive GSTC-FDTD results match well with the analytical values, with no respect to the FDTD grid size.
As the next example, a two-dimensional problem for a BP sheet is analyzed and the geometry of this problem is shown in Fig. 7. We consider two different BP sheet orientations. For  the armchair case, we set the electron mass having armchair characteristic along the x direction (and zigzag characteristic along the z direction). For the zigzag case, we set the electron mass having zigzag characteristic along the x direction (and VOLUME 10, 2022  armchair characteristic along the z direction). We employ the FDTD grid size of 1 nm, 2 nm, and 10 nm and the corresponding time step size is set to 2.333×10 −18 s, 4.667×10 −18 s, and 2.333 × 10 −17 s respectively for the CFL stability condition. The physical size of the FDTD domain is 1000 nm × 1000 nm and 10-cell PML layers are used in the each direction. The length of the BP sheet along the x direction is 1000 nm. The z-directed magnetic current with the differentiated Gaussian pulse in the time domain is excited at 20 nm above the BP  sheet and then H z at 10 nm below the BP sheet is observed. Other simulation setup parameters are maintained same as in the one-dimensional case. Figs. 8, 9, and 10 shows the normalized magnetic fields of the three numerical results  when the FDTD grid size is 1 nm, 2 nm, and 10 nm. The numerical result of both dispersive GSTC-FDTDs are same, regardless of the FDTD grid size. However, the conventional FDTD results are in agreement with dispersive GSTC-FDTD results for only the FDTD grid size of 1 nm (Fig. 8) and thus the conventional FDTD simulations with larger FDTD sizes (Figs. 9 and 10) cannot be used. Table 2 summarizes CPU time and memory storage for accurate electromagnetic simulations: both GSTC-FDTDs with the FDTD grid size of 10 nm and the conventional FDTD with the FDTD grid size of 1 nm. As shown in the table, the proposed Drudedipsersive GSTC-FDTD simulation requires less CPU time and memory storage requirement compared to the other two FDTD simulations.

V. CONCLUSION
In this work, we propose an efficient dispersive GSTC-FDTD formulation based on the Drude dispersion model. The number of field variables and operation count is discussed for the proposed dispersive GSTC-FDTD formulation and the previous Lorent-based GSTC-FDTD formulation. In addition, the numerical surface susceptibilities of the proposed Drude-dispersive GSTC-FDTD algorithm and the conventional Lorentz-dispersive GSTC-FDTD algorithm are derived to confirm the numerical accuracy. Numerical examples are employed to demonstrate that the proposed Drude-dispersive GSTC-FDTD formulation can accurately and effectively analyze electromagnetic waves interaction with the metasurface.