Prismatic discontinuous Galerkin time domain method with an integrated generalized dispersion model for efficient optical metasurface analysis

Wending Mai; Sawyer D. Campbell; Eric B. Whiting; Lei Kang; Pingjuan L. Werner; Yifan Chen; Douglas H. Werner

doi:10.1364/OME.399414

1. Introduction

Optical metamaterials are a type of engineered electromagnetic material that interacts with light at terahertz (THz), infrared (IR), or visible wavelengths. Optical metasurfaces comprise a class of optical metamaterials with subwavelength thicknesses and exceptional abilities for controlling light. Such two-dimensional (2D) and quasi-2D metamaterials provide designers with the capability to fully control light with planar elements and opens the door to a new generation of planar photonics devices [1]. An important milestone in planar photonics is the recent development metasurface design techniques which have brought transformative functionality and enhanced performance with reduced dimensions compared to conventional optical elements. As such, many metasurface designs are compatible with existing planar, low-cost manufacturing technologies in modern industry. Furthermore, some metasurfaces are also compatible with planar electronic devices, which is of critical importance for applications in opto-electronics, microscopy, imaging, and sensing [2,3]. In order to exploit the advantages of planar photonics, a sophisticated designed algorithm with high accuracy and efficiency is favored for numerical simulation.

Time-domain computational electromagnetic solvers have advantages for providing transient information and wideband responses. The Finite Difference Time Domain (FDTD) method is the most well-known transient solver for its high efficiency and simplicity. Comparatively, the Discontinuous Galerkin Time Domain (DGTD) inherently possesses both superior meshing (h-) and high-order (p-) accuracy [4–8], which are intrinsic disadvantages of FDTD-based solvers. Moreover, the DGTD method can be regarded as an element-level decomposition technique; all operations of DGTD are local and easily parallelizable. It has become a popular choice for solving transient problems as it integrates the advantages of finite volume time domain method (FVTD) [9] and finite element time domain method (FETD) [10] in an explicit time-marching scheme. In the terahertz, infrared, and optical frequency regimes, material dispersion in metals cannot be neglected. In order the accurately model the frequency dependent permittivity of metals in time domain methods, researchers have developed various fitting models. The Drude-critical point (DCP) model has shown its own superiority over Drude and Drude-Lorentz models to approximate the dispersion of metals. Further still, a generalized dispersive material (GDM) model based on [1/2] Padé approximants, which covers Debye [11–15], Drude [13,16–18], Lorentz [13], Sellmeier [19], and critical points models [20] with an arbitrary number of poles, has been proposed as a robust universal model for materials with arbitrary dispersive behaviors [16]. Recently, the GDM model was incorporated into DGTD to facilitate full wave time domain simulation of optical devices that contain dispersive materials [21]. Furthermore, periodic boundary conditions (PBCs) were implemented for analysis of periodic nanostructures under oblique incidence [22].

In practice, metasurfaces are usually represented by a patterned metallo-dielectric layer that is very thin compared to the wavelength of the incident light [23–26]. For low-cost manufacturing, most metasurfaces are composed of planar structures. The treatment of this thin layer is challenging for conventional numerical methods as it leads to a multi-scale modeling problem [27]. One option is to treat the metal layer as a type of boundary condition without thickness. This includes surface impedance boundary conditions (SIBC) and the generalized sheet transition conditions (GSTC) [23,24]. However, neglecting the layer thickness will bring unexpected approximation error. Proper meshing of the thin layer is required for better accuracy. Tetrahedral mesh elements are the most popular but may be a poor choice and lead to low quality when discretizing such thin structures. Of course, low quality meshes will bring poor convergence of the algorithm, or even inaccurate results. While hexahedral mesh elements do not suffer from this issue, they may be lacking in their ability to accurately describe complex geometrics. On the contrary, triangular prism meshing is advantageous when compared with other conventional meshing techniques for planar metasurfaces as it can accurately model arbitrary patterns while being free of the quality issue that limits tetrahedral elements, even for ultra-thin layered geometries. The prismatic DGTD method has been utilized to analyze high-speed printed circuit boards (PCB) and microwave frequency selective surfaces (FSS) [28–30]. However, there are currently no available DGTD schemes that are capable of efficiently analyzing thin dispersive metallo-dielectric optical metasurfaces.

In this paper, the GDM model is integrated with prism-based DGTD for efficient simulation of planar photonic devices, e.g. metasurfaces, comprised of dispersive metallo-dielectric structures in the optical regime. The accuracy and convergence of the proposed prism-based DGTD + GDM method is validated by comparing it with analytical results. A metasurface composed of a thin gold film with a periodic array of air holes, and a pixelized gold metasurface were analyzed to validate the accuracy and efficiency of the proposed algorithm. In Section 2, the prism-based DGTD + GDM formulation is presented. Several numerical examples are provided in Section 3 to validate the accuracy and efficiency improvement of the proposed method.

2. Methods and formulation

2.1 Governing equations for the generalized dispersive model

For a linear isotropic homogeneous medium, the dispersive form of Maxwell’s symmetric curl equations are given by:

(1)$${\mu _0}{\mu _r}j\omega \vec{{\boldsymbol H}} ={-} \nabla \times \vec{{\boldsymbol E}}$$

(2)$${\varepsilon _0}{\varepsilon _r}j\omega \vec{{\boldsymbol E}} = \nabla \times \vec{{\boldsymbol H}}$$

where $\vec{{\boldsymbol E}}$ and $\vec{{\boldsymbol H}}$ represent the frequency domain (bold) electric and magnetic field vectors, respectively, while the component $j\omega $ in the frequency domain represents the engineering convention corresponding to time-varying fields as ${e^{j\omega t}}$.

In order to fit the dispersive relative permittivity in time domain, here we employ the GDM model which can be expressed as a sum of the permittivity ${\varepsilon _r}$ at the infinite high frequency limit (${\varepsilon _\infty })$ and m partial terms [20]:

(3)$${\varepsilon _r} = {\varepsilon _\infty } + \mathop \sum \limits_m \frac{{{a_{0,m}} + j\omega {a_{1,m}}}}{{{b_{0,m}} + j\omega {b_{1,m}} + {\omega ^2}}}$$

The performance of the GDM model is dependent on its associated frequency independent parameters ${a_0}$, ${a_1}$, ${b_0}$, and ${b_1}$ (${b_0} \ge 0$). The GDM model can be rearranged to represent multiple classical dispersion formulas such as the Drude (${a_{1,m}} = {b_{0,m}} = 0)$, Lorentz (${a_{1,m}} = 0)$, Sellmeier (${a_{1,m}} = {b_{1,m}} = 0)$, and critical points models [20,31].

The generalized dispersion model (GDM) comes from the Padé approximation of an argument $- j\omega $ with real coefficients. Therefore, it does not naturally conform to causal physics with general parameters. Thus, it does not always satisfy the Kramers-Kronig (KK) relations for arbitrary parameters. However, in many applications, the parameters of the GDM can be selected to conform with causality. For example, in this paper, we used the D2CP model which contains a Drude term and two critical point terms. For the critical point model, the conventional approach is to first fit the imaginary part, then a match of the real parts automatically results from KK-consistency. For the Drude term containing real poles, the frequency domain dielectric function must have additional delta-function terms to comply with the KK relation. However, it is convention that all real poles are far offset from the intended frequency range, and are therefore omitted. To summerize, the GDM does not essentially satisfy the KK relation, but the D2CP model we used in this paper and its parameter settings does satisfy the KK relation in the frequency range of interest [31].

By substituting (3) into (2), we have:

(4)$$\nabla \times \vec{{\boldsymbol H}} = j\omega ({\varepsilon _0}{\varepsilon _\infty }\vec{{\boldsymbol E}} + {\varepsilon _0}\mathop \sum \limits_m {\overrightarrow {\boldsymbol P} _m})$$

where the introduced polarization vector ${\overrightarrow {\boldsymbol P} _m}$ is defined as:

(5)$${\overrightarrow {\boldsymbol P} _m} = \frac{{{a_{0,m}} + j\omega {a_{1,m}}}}{{{b_{0,m}} + j\omega {b_{1,m}} - {\omega ^2}}}\vec{{\boldsymbol E}}$$

Through application of an inverse Fourier transformation, (5) can be recast in the form of an auxiliary differential equation (ADE) [13,16,21,22]:

(6)$$\frac{{{d^2}{{\overrightarrow {\boldsymbol P} }_m}}}{{d{t^2}}} + {b_{1,m}}\frac{{d{{\overrightarrow {\boldsymbol P} }_m}}}{{dt}} + {b_{0,m}}{\overrightarrow {\boldsymbol P} _m} = {a_{0,m}}\vec{E} + {a_{1,m}}\frac{{d\vec{E}}}{{dt}}$$

Accordingly, the dispersive form of Maxwell’s equations in (1) and (2) can be transformed into the time domain as:

(7)$${\mathrm{\mu}_0}{\mathrm{\mu}_r}\frac{{d\vec{H}}}{{dt}} ={-} \nabla \times \vec{E}$$

(8)$${{\varepsilon }_0}{\varepsilon }_\infty \frac{{d\vec{E}}}{{dt}} + {{\varepsilon }_0}\mathop \sum \limits_m \frac{{d{{\overrightarrow {\boldsymbol P} }_m}}}{{dt}} = \nabla \times \vec{H}$$

where $\vec{E}$, $\vec{H}$, and ${\overrightarrow {\boldsymbol P} _m}$ are the electric, magnetic and polarization vectors in the time domain (unbold).

2.2 Spatial discretization of the prism DGTD

Fig. 1. Geometrical illustration of a prismatic element and its associated basis functions.

Download Full Size | PDF

In (7), the electric and magnetic fields, as well as the electric polarization vector are expanded by the basis functions $\vec{\Phi }_k^i({\vec{r}} )$ in every prismatic element $ {\Omega _i}$:$\overrightarrow {{E^i}} ({\vec{r},t} )= \mathop \sum \nolimits_{k = 1}^{n_e^i} e_k^i(t )\vec{\Phi }_k^i({\vec{r}} )$, $\overrightarrow {{H^i}} ({\vec{r},t} )= \mathop \sum \nolimits_{l = 1}^{n_h^i} h_l^i(t )\vec{\Phi }_l^i({\vec{r}} )$, and $\overrightarrow {P_m^n} ({\vec{r},t} )= \mathop \sum \nolimits_{n = 1}^{n_p^i} p_{n,m}^i(t )\vec{\Phi }_n^i({\vec{r}} )$ with $n_e^i$, $n_h^i$, and $n_p^i$ representing the number of basis functions for $\vec{E}$, $\vec{H}$ and $\vec{J}$ in element i, respectively.

Accordingly, the vector basis functions $\vec{\Phi }_k^i({\vec{r}} )$ in Fig. 1 are are defined on the triangular prism element as [32,33]:

(9)$$\left\{ \begin{array}{l} {\overrightarrow \Phi _\textrm{1}}\textrm{ = }\zeta ({{L_\textrm{1}}\nabla {L_\textrm{2}}\textrm{ - }{L_\textrm{2}}\nabla {L_\textrm{1}}} )\\ {\overrightarrow \Phi _\textrm{2}}\textrm{ = }\zeta ({{L_\textrm{2}}\nabla {L_\textrm{3}}\textrm{ - }{L_\textrm{3}}\nabla {L_\textrm{2}}} )\\ {\overrightarrow \Phi _\textrm{3}}\textrm{ = }\zeta ({{L_\textrm{3}}\nabla {L_\textrm{1}}\textrm{ - }{L_\textrm{1}}\nabla {L_\textrm{3}}} )\\ {\overrightarrow \Phi _\textrm{4}}\textrm{ = }({\textrm{1 - }\zeta } )({{L_\textrm{4}}\nabla {L_\textrm{5}}\textrm{ - }{L_\textrm{5}}\nabla {L_\textrm{4}}} )\\ {\overrightarrow \Phi _\textrm{5}}\textrm{ = }({\textrm{1 - }\zeta } )({{L_\textrm{5}}\nabla {L_\textrm{6}}\textrm{ - }{L_\textrm{6}}\nabla {L_\textrm{5}}} )\\ {\overrightarrow \Phi _\textrm{6}}\textrm{ = }({\textrm{1 - }\zeta } )({{L_\textrm{6}}\nabla {L_\textrm{4}}\textrm{ - }{L_\textrm{4}}\nabla {L_\textrm{6}}} )\\ {\overrightarrow \Phi _\textrm{7}}\textrm{ = }{L_\textrm{1}}\nabla \zeta \textrm{/|}\nabla \zeta \textrm{|}\\ {\overrightarrow \Phi _\textrm{8}}\textrm{ = }{L_\textrm{2}}\nabla \zeta \textrm{/|}\nabla \zeta \textrm{|}\\ {\overrightarrow \Phi _\textrm{9}}\textrm{ = }{L_\textrm{3}}\nabla \zeta \textrm{/|}\nabla \zeta \textrm{|} \end{array} \right.$$

where ${L_i}$ denotes the nodal basis function at node i with i=1, …, 6:

(10)$${L_i} = \frac{1}{{6{V^e}}}({a_i^e + b_i^ex + c_i^ey + d_i^ez} )$$

The term ${V^e}$ represents the volume of the element, the coefficients $a_i^e,b_i^e,c_i^e$ and $d_i^e$ can be determined the same way as in [10]. The parameter $\zeta $ varies from 0 to 1 along the height of the prism. It is noted that the basis functions in the vertical direction ${\vec{\Phi }_i}$, i = 7, …, 9 are orthogonal to the basis functions in the horizontal direction ${\vec{\Phi }_i}$, i = 1, . . ., 6 [12].

To model the curl equation in (7) and (8), we introduce a numerical upwind flux that is based on the Rankine Hugoniot jump relations. The numerical flux is implemented via the weak boundary condition of DGTD [5,8]:

(11)$$\hat{n}_i^f \times \vec{H}_f^{\ast } = \hat{n}_i^f \times \frac{{({{Z^i}{{\vec{H}}^i} + Z_f^j\vec{H}_f^j} )+ \hat{n}_i^f \times ({{{\vec{E}}^i} - \vec{E}_f^j} )}}{{{Z^i} + Z_f^i}} - \beta \frac{{\hat{n}_i^f \times {{\vec{M}}_s}}}{{({{Z^i} + Z_f^j} )}}$$

(12)$$\hat{n}_i^f \times \vec{E}_f^{\ast } = \hat{n}_i^f \times \frac{{({{Y^i}{{\vec{E}}^i} + Y_f^j\vec{E}_f^j} )+ \hat{n}_i^f \times ({\vec{H}_f^j - {{\vec{H}}^i}} )}}{{{Y^i} + Y_f^i}} - \beta \frac{{\textrm{Y}_f^j{{\vec{M}}_s}}}{{({{Y^i} + Y_f^j} )}}$$

where j represents the neighboring element adjacent to the face f of element i, $\hat{n}_i^f$ denotes the outward unit vector to face f, $\hat{n}_i^f \times \vec{H}_f^{\ast }$ and $\hat{n}_i^f \times \vec{E}_f^{\ast }$ are the numerical fluxes used for the exchange between adjacent elements required by the weak boundary condition. ${Z^{i/j}} = \sqrt {{\mu ^{i/j}}/{\varepsilon ^{i/j}}} $ and ${Y^{i/j}} = 1/{Z^{i/j}}$ denote the characteristic impedance and admittance, respectively. Finally, we note that $\beta = 1$ if the face f is over the wave port, otherwise $\beta = 0$.

It should be pointed out that the impedances ${Z^{i/j}}$ in the upwind flux (11) and (12) contain the material permittivity ${\varepsilon }(\mathrm{\omega} )$, and thus are inherently dispersive. Consequently, the solution of (11) and (12) should be convolutional. However, in practice, we could choose ${Z^{i/j}} = \sqrt {{\mu _\infty }^{i/j}/{\varepsilon _\infty }^{i/j}} $ as a sufficient approximation to avoid the need to carry out a time-consuming convolution [24].

By performing Galerkin testing over (7) and (8), and taking into account the numerical flux in (11) and (12), we obtain the following DGTD semi-discretized matrix equations:

(13)$$\begin{array}{l} \frac{{d{e^i}}}{{dt}} = \{ \bar{M}{_e^{i - 1}}\left[ {\bar{S}_e^i{h^i} + \mathop \sum \limits_{f = 1}^{n_i^f} ({\bar{F}_{ee}^{ii,f}e_f^i + \bar{F}_{ee}^{ij,f}e_f^j + \bar{F}_{eh}^{ii,f}h_f^i + \bar{F}_{eh}^{ij,f}h_f^j} )+ \beta \bar{F}_e^{i,{M_s}}} \right]\\ - {\varepsilon _0}\mathop \sum \limits_m \frac{{\mathrm{\partial} p_m^i}}{{\mathrm{\partial} t}}\} /({{\varepsilon_0}{\varepsilon_\infty }} )\end{array}$$

(14)$$\frac{{d{h^i}}}{{dt}} = \bar{M}{_h^{i - 1}}\left[ { - \bar{S}_h^i{e^i} + \mathop \sum \limits_{f = 1}^{n_i^f} ({\bar{F}_{hh}^{ii,f}h_f^i + \bar{F}_{hh}^{ij,f}h_f^j + \bar{F}_{he}^{ii,f}e_f^i + \bar{F}_{he}^{ij,f}e_f^j} )+ \beta \bar{F}_h^{i,{M_s}}} \right]/{\mu _0}$$

where $\bar{M},\bar{S},$ and $\bar{F}$ denote the mass matrix, the stiffness matrix and the flux matrix, respectively, whose detailed definitions are given as:

{[{\bar{M}_e^i} ]_{kl}} = \mathop \smallint \nolimits_{{\Omega _i}}^{} \vec{\Phi }_k^i \cdot {\varepsilon ^i}\vec{\Phi }_l^id\vec{r}

{[{\bar{M}_h^i} ]_{kl}} = \mathop \smallint \nolimits_{{\Omega _i}}^{} \vec{\Phi }_k^i \cdot {\mu ^i}\vec{\Phi }_l^id\vec{r}

{[{\bar{M}_q^i} ]_{kl}} = {[{\bar{M}_{c,q}^i} ]_{kl}} = \mathop \smallint \nolimits_{{\Omega _i}}^{} \vec{\Phi }_{k,q}^i \cdot \vec{\Phi }_{l,q}^id\vec{r}

{[{\bar{S}_e^i} ]_{kl}} = {[{\bar{S}_h^i} ]_{kl}} = \mathop \smallint \nolimits_{{\Omega _i}}^{} \vec{\Phi }_k^i \cdot \nabla \times \vec{\Phi }_l^id\vec{r}

{[{\bar{F}_{ee}^{ii,f}} ]_{kl}} = \frac{{ - 1}}{{{Z^i} + Z_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot \hat{n}_i^f \times ({\hat{n}_i^f \times \vec{\Phi }_l^i} )d\vec{r}

{[{\bar{F}_{ee}^{ij,f}} ]_{kl}} = \frac{1}{{{Z^i} + Z_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot \hat{n}_i^f \times ({\hat{n}_i^f \times \vec{\Phi }_l^{j,f}} )d\vec{r}

{[{\bar{F}_{eh}^{ii,f}} ]_{kl}} = \frac{{ - Z_f^i}}{{{Z^i} + Z_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot \hat{n}_i^f \times \vec{\Phi }_l^id\vec{r}

{[{\bar{F}_{eh}^{ij,f}} ]_{kl}} = \frac{{Z_f^i}}{{{Z^i} + Z_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot \hat{n}_i^f \times \vec{\Phi }_l^{j,f}d\vec{r}

{[{\bar{F}_e^{i,{M_s}}} ]_k} = \frac{{ - 1}}{{{Z^i} + Z_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot ({\hat{n}_i^f \times {{\vec{M}}_s}} )d\vec{r}

{[{\bar{F}_{hh}^{ii,f}} ]_{kl}} = \frac{{ - 1}}{{{Y^i} + Y_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot \hat{n}_i^f \times ({\hat{n}_i^f \times \vec{\Phi }_l^i} )d\vec{r}

{[{\bar{F}_{hh}^{ij,f}} ]_{kl}} = \frac{1}{{{Y^i} + Y_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot \hat{n}_i^f \times ({\hat{n}_i^f \times \vec{\Phi }_l^{j,f}} )d\vec{r}

{[{\bar{F}_{he}^{ii,f}} ]_{kl}} = \frac{{Y_f^i}}{{{Y^i} + Y_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot \hat{n}_i^f \times \vec{\Phi }_l^id\vec{r}

{[{\bar{F}_{he}^{ij,f}} ]_{kl}} = \frac{{ - Y_f^i}}{{{Y^i} + Y_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot \hat{n}_i^f \times \vec{\Phi }_l^{j,f}d\vec{r}

(15)$${[{\bar{F}_h^{i,{M_s}}} ]_k} = \frac{{ - Y_f^i}}{{{Y^i} + Y_f^j}}\mathop \smallint \nolimits_{\mathrm{\partial} \Omega _i^f}^{} \vec{\Phi }_k^i \cdot {\vec{M}_s}d\vec{r}$$

Also, ${\bar{M}^{ - 1}}$ in (13) and (14) is the inverse matrix corresponding to $\bar{M}$. The quantities ${e^i}$ and ${h^i}$ are the electronic and magnetic column vectors containing the unknown coefficients in element i. The coefficients with a superscript f correspond to those on face f between element i and j.

2.3 Runge-Kutta time stepping scheme

Next, we adopt a time integration scheme to solve (6), (13) and (14) with a temporal discretization that meets the Courant-Freidrichs-Lewy (CFL) constraint. The low-storage explicit Runge-Kutta method (LSERK4) is widely used in DGTD for its lower memory consumption compared to the fourth-order, four-stage explicit Runge-Kutta method (ERK) (ERK4). However, studies have shown that for dispersive materials, the introduced auxiliary differential function (ADE) will impact the eigen value of the governing equations, and eventually cause a stability issue with LSERK4 scheme [34]. Based on this observation, most of the current works on dispersive materials use the ERK4 or the LF2 scheme for stability considerations [13,22]. In this paper, we choose the ERK4 scheme for time integration.

In order to apply the ERK4 method to solve the second order derivative term in (6), an auxiliary vector is introduced to lower the derivative order, which is also decomposed into the prismatic basis functions defined in (9), such that $\overrightarrow {Q_m^n} ({r,t} )= \mathop \sum \nolimits_{n = 1}^{n_p^i} q_{n,m}^i(t )\vec{\Phi }_n^i(r )$. Thus, we have:

(16)$$q_m^i = \frac{{dp_m^i}}{{dt}}$$

In such, the Eqs. (13) and (6) can be rewritten as:

(17)$$\begin{array}{l} \frac{{d{e^i}}}{{dt}} = \{ \bar{M}{_e^{i - 1}}\left[ {\bar{S}_e^i{h^i} + \mathop \sum \limits_{f = 1}^{n_i^f} ({\bar{F}_{ee}^{ii,f}e_f^i + \bar{F}_{ee}^{ij,f}e_f^j + \bar{F}_{eh}^{ii,f}h_f^i + \bar{F}_{eh}^{ij,f}h_f^j} )+ \beta \bar{F}_e^{i,{M_s}}} \right]\\ - {\varepsilon _0}\mathop \sum \limits_m q_m^i\} /({{\varepsilon_0}{\varepsilon_\infty }} )\end{array}$$

(18)$$\frac{{dq_m^i}}{{dt}} + {b_{1,m}}q_m^i + {b_{0,m}}p_m^i = {a_{0,m}}{e^i} + {a_{1,m}}\frac{{d{e^i}}}{{dt}}$$

Then, the fourth-order four-stage explicit Runge-Kutta method (ERK) is adopted by setting the operator ${{{\cal L}}_i}({u_n^i,{t_n}} )$ equal to the right-hand side of (14), (17) and (18). These three first-order differential equations can be evaluated by using the following expressions:

(19)$$\left\{ \begin{array}{l} {k^{(1)}} = {{{\cal L}}_i}({u_n^i,{t_n}} )\\ {k^{(2)}} = {{{\cal L}}_i}\left( {u_n^i + \frac{1}{2}\delta t\ast {k^{(1)}},{t_n} + \frac{1}{2}\delta t} \right)\\ {k^{(3)}} = {{{\cal L}}_i}\left( {u_n^i + \frac{1}{2}\delta t\ast {k^{(2)}},{t_n} + \frac{1}{2}\delta t} \right)\\ {k^{(4)}} = {{{\cal L}}_i}({u_n^i + \delta t\ast {k^{(3)}},{t_n} + \delta t} )\\ u_{n + 1}^i\textrm{ = }u_n^i + \frac{1}{6}\delta t\ast ({{k^{(1)}} + 2{k^{(2)}} + 2{k^{(3)}} + {k^{(4)}}} )\end{array} \right.$$

where $u_n^i$ represents the unknowns $e_n^i$, $h_n^i$, $p_n^i$ or $q_n^i$ when solving for the electric field, magnetic field, polarized current vector and the auxiliary vector, respectively. Also, $\delta t$ is the maximum time-step size for the DGTD scheme, which is determined by the CFL condition.

2.4 Extraction of transmission and reflection

To make use of periodic metasurface structures with ${D_{2n}}$ symmetry, we assigned a pair of PEC and PMC boundaries to reduce the number of unknowns. In this way, only a quarter of the structure is computed and thus it is more efficient. Accordingly, the normally incident plane wave is excited by a TEM mode waveguide port whose field distribution is exactly the same as a plane wave [30].

Transmission and reflection are the typically used metrics to demonstrate the performance of an optical metasurface. To extract the transmission and reflection data, we first expand the electric coefficient corresponding to the TEM mode at the input and output wave port.

(20)$$a_{TEM}^j = \smallint \vec{E} \cdot \overrightarrow {{e_{TEM}}} dS$$

where j denotes the input (i) and output port (o), respectively. Moreover, $\overrightarrow {{e_{TEM}}} $ is the TEM mode distribution on the wave port. Then the reflection (R) and transmission (T) data over the given wave ports can be conveniently evaluated by:

(21)$$\left\{ {\begin{array}{l} {R(t )= \frac{{a_{TEM}^i - \tilde{a}_{TEM}^i}}{{\tilde{a}_{TEM}^i}}\; \; \; \; }\\ {T(t )= \frac{{a_{TEM}^o}}{{\tilde{a}_{TEM}^i}}\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; } \end{array}} \right.$$

where $\tilde{a}_{TEM}^j$ is the mode expansion coefficient of the incident electric field. Finally, a Fourier transformation can be used to convert the transmission and reflection from the time domain (21) to the frequency domain:

(22)$$\left\{ {\begin{array}{l} {|R |= \left|{FFT\left( {\frac{{a_{TEM}^i - \tilde{a}_{TEM}^i}}{{\tilde{a}_{TEM}^i}}} \right)} \right|\; \; \; \; }\\ {|T |= \left|{FFT\left( {\frac{{a_{TEM}^o}}{{\tilde{a}_{TEM}^i}}} \right)} \right|\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; } \end{array}} \right.$$

3. Numerical examples

In this section, various benchmarks are given to validate the accuracy, flexibility, and efficiency of the proposed algorithm. All examples employ a magnetic current excitation with amplitude of ${A_0}(t )$, described by a sinusoidally modulated Gaussian pulse defined as:

(23)$${A_0}(t )={-} \exp \left[ { - \frac{{4\mathrm{\pi}{{({t - {t_0}} )}^2}}}{{{t_n}}}} \right]\textrm{cos}[{2\mathrm{\pi}ft} ]$$

where ${t_0}$ denotes the reference temporal points, ${t_n}$ is a spectral bandwidth related parameter, and f is the frequency of the modulated signal.

All numerical simulations were performed on a laptop with an Intel i6700HQ CPU of 2.6 GHz, 4 cores, 8 threads, and 16 GB of memory. The algorithm has been fully parallelized with all 8 threads by OpenMP.

3.1 Optimal settings of the generalized dispersive model

The optical properties of gold are difficult to represent in the visible/near-UV region with an analytic model. The reason for this is the more important role played by interband transitions in gold in the violet/near-UV region. Gold has at least two interband transitions at about 470 and 330 nm that do play an important role and must be included explicitly for a realistic analytic model [35]. Consequently, for λ > 500 nm, a Drude model is sufficient for representing the dispersion of gold. However, a more sophisticated analytic model has to be employed for smaller wavelengths.

To accurately describe the dispersion of gold over a broad bandwidth, here we employed a specific kind of GDM composed of one Drude term and two critical point terms (D2CP) [11,29]. As mentioned, the GDM model can be represented as a summation of classical dispersion models. The relationship of GDM parameters ${a_0}$, ${a_1}$, ${b_0}$, and ${b_1}$ with other dispersion model parameters can be found in [20]. As a result, the parameters of the D2CP model for gold can be readily derived and are provided in Table 1.

Table 1. Parameters of the GDM model composed of D2CP

View Table | View all tables in this article

Figure 2(a) demonstrates the real and imaginary parts of the relative permittivity of the D2CP GDM model. It can be seen that the D2CP GDM yields a good fit to the experimental measurements over an ultra-wide band [36]. As such, the proposed prism-based DGTD formulation with an integrated GDM model is capable of broadband simulation for dispersive metals in the optical regime.

Fig. 2. Comparison between the relative permittivity of gold from the D2CP GDM model and experimentally determined results. The squares represent the real part of the experimental values, while the diamonds correspond to the imaginary part [33]. The solid and dashed lines represent the real and imaginary parts of the corresponding D2CP GDM model, respectively. The D2CP GDM model exhibits good accuracy over a broad bandwidth.

Download Full Size | PDF

3.2 Reflection and transmission of a thin gold film

We next study the reflection and transmission of a thin gold film upon normal incidence. The analytical result can be derived from closed-formed Fresnel equations. By comparing the numerical results of the proposed algorithm with the analytical solution, we can determine the accuracy and convergence rate of the proposed method.

The gold film has a thickness of 5 nm, which is ultra-thin in contrast with the minimum assumed wavelength of 188 nm. The film is illuminated by a sinusoidally modulated Gaussian pulse defined in (23). The transient response of the reflection and transmission coefficients is shown in Fig. 3(a) with a mesh configuration consisting of prismatic elements with a maxima length of 12 nm. Through Fourier transformation, both the reflection and transmission coefficients can be recovered in the chosen wavelength range. As shown in Fig. 3(b), the result of the proposed DGTD + GDM algorithm matches very well with the analytical data. As can be seen, a thin gold film is almost transparent in the visible spectrum, but becomes more reflective as the wavelength increases, e.g. for the infrared and terahertz regimes.

Fig. 3. The reflection and transmission properties of a thin gold film. (a) The transient response. (b) The reflection and transmission results from the proposed prism DGTD + GDM method and from the analytical expression. (c) Convergence plot showing accuracy improvement along with the refined mesh size.

Download Full Size | PDF

In order to investigate the accuracy of the proposed DGTD + GDM method, the relative error of the numerical and analytical results are depicted as a function of the mesh size. Figure 3(c) shows the convergence plot which provides a means for quantifying the degree of accuracy improvement against refinements in the mesh [4,14–18]. The algorithm was tested for different mesh sizes, where a measure of the accuracy was determined from the root mean square error (RMSE), which here is defined as:

(24)$$RMSE = \sqrt {\frac{{{{\sum\nolimits_1^{{N_\lambda }} {(R/T} }^{DGTD + GDM}} - R/{T^{Analytical}}{)^2}}}{{{N_\lambda }}}}$$

such that the performance was compared at ${N_\lambda }$ sampling wavelengths.

3.3 Illumination on a thin gold hole array

At this point we apply the prism-based DGTD method with the GDM model outlined in the previous section to compute the S-parameters of a representative thin gold hole array structure for further validation. Figure 4 illustrates a single unit cell of the doubly periodic infinite array. The unit cell of the hole array under consideration contains a thin gold film layer with a thickness of 30 nm. The 100 nm radius holes are spaced by a 400 nm lattice period.

Fig. 4. Comparison of the gold hole array unit cell spatially discretized as (a) prisms and (b) tetrahedrons.

Download Full Size | PDF

Figure 4(a) also depicts the prismatic spatial discretization utilized by the proposed algorithm, while Fig. 4(b) represents the conventional tetrahedral spatial discretization. As can be seen, with similar element lengths, the prismatic mesh cells allow for fewer elements, which bring with them fewer unknowns and higher computational efficiency.

Figure 5 depicts the computational domain of the numerical example. To make use of the periodic metasurface structures with ${D_{2n}}$ symmetry, we assigned a pair of PEC and PMC boundaries in the x- and y-directions, respectively. In this way, only a quarter of the structure is computed, which improved the computational efficiency. The normally incident plane wave is excited by a TEM mode waveguide port at -500 nm. The computational domain is truncated by a pair of Silver-Mueller absorbing (SMA) boundary conditions in the z- direction at -700 nm and 700 nm, respectively. The SMA boundary condition provides an ideally null reflection coefficient for normal incidence [34].

Fig. 5. The computational domain of the prism DGTD method. It is truncated by a pair of PEC, PMC and SBC boundary conditions in x-, y- and z- directions, respectively.

Download Full Size | PDF

The minimum element length is about 20 nm, while the time step δt is set to 4.2 atto seconds. The transient results generated by the prism-based DGTD with the GDM model are shown in Fig. 6(a). A total of 90,000 time steps were computed, corresponding to nearly 37 femto seconds. In order to demonstrate the frequency selective properties of the nanohole array, the transient S-parameters shown in Fig. 6(a) are Fourier transformed to the frequency domain and plotted in Fig. 6(b). For comparison, Fig. 6(b) also contains the simulated result obtained from the conventional tetrahedral mesh and from the commercial CST software package. The comparison was made with simulation results from CST’s frequency domain solver (FEM) using imported experimental dispersion data computed at 20 frequency points [36]. For a fair comparison, the PEC + PMC boundaries setting was also adopted in the CST solver. Figure 6(b) shows that the result of the proposed algorithm yields good agreement with the CST FEM solver. Some difference exists because of the discrepancy between the experimental material dispersion and the D2CP GDM model.

Fig. 6. (a) The transient response of a gold nano-hole array. (b) Fourier transformed results in the wavelength regime compared with results of a conventional tetrahedral mesh and commercial software.

Download Full Size | PDF

Moreover, the prism-based DGTD with the GDM algorithm requires only 230 seconds to perform the computations, while the CST FEM solver consumes 2,658 seconds. This computational efficiency enhancement is a direct consequence of the optimal spatial discretization enabled by prismatic elements, and the ease by which parallel computing can be utilized within the DGTD framework.

The gold nano-hole array demonstrates a remarkable resonance observed at 680 nm. As presented in [37], this type of extraordinary transmission behavior can be primarily attributed to the excitation of a surface plasmon at the metallic hole array structure interface. Concisely, the incident light couples into electromagnetic surface waves (i.e., surface plasmon polaritons (SPPs) at the metal-dielectric interface), which then radiate through reciprocal interactions with the structure. This, in turn, produces unique features in the transmission and reflection spectra.

The most advantageous aspect of the proposed prism DGTD + GDM algorithm is its efficiency improvement compared to other solvers. In Table 2, we compared the number of mesh cells, the degree of freedom, and the computation CPU time for different methods. It can be seen that the prismatic meshing technique has fewer mesh cells than the tetrahedral mesh, given the same element length. Consequently, it has fewer unknowns and higher computational efficiency. As shown, the DGTD methods with both tetra and prism meshes have higher efficiency compared to the CST FEM solver because of the transient computation and element level parallelization. Moreover, the prismatic mesh has fewer elements than the tetrahedral mesh, given the same element length. The tetrahedral and the triangular prism mesh elements have 6 and 9 unknowns for every component, respectively. Except for the electric and magnetic field ($\vec{E}$ and $\vec{H}$), the auxiliary components of ${\vec{P}_{1,2,3}}$ and ${\vec{Q}_{1,2,3}}$ are also computed. In total, the prism DGTD + GDM method has significantly fewer unknowns than a conventional tetrahedral DGTD + GDM method. Consequently, it achieves significantly higher computational efficiency, namely, it 4 times faster than conventional methods.

Table 2. Comparison of the number of elements, number of unknowns, CPU time, and memory consumption for different methods^a

View Table | View all tables in this article

We also notice a reduction in the memory consumption for the prism DGTD + GDM when compared with the tetra DGTD + GDM. This memory reduction comes from an optimal space discretization which requires fewer elements and fewer unknowns. However, both the tetra and prism DGTD methods have higher memory consumption than the CST FEM solver. This partly originates from the extra storage requirement for the auxiliary components (${\vec{P}_{1,2,3}}$ and ${\vec{Q}_{1,2,3}}$) and the overheads of OpenMP parallel computing.

3.4 Pixelized optical metasurface

Among various metasurface designs, pixelized metasurfaces represent one of the most popular classes as they possess a high degree of design freedom [38,39]. In order to exploit the advantages of these pixelized metasurfaces, many advanced optimization algorithms have been developed to ensure that a successful design can be found with the desired functionality [40]. An accurate and highly efficient full-wave simulation technique is an essential part of the optimization process, and is especially important for pixelized metasurface design.

We apply the proposed algorithm to analyze a representative pixelized metasurface which is composed of a doubly-infinite array of unit cells (see Fig. 7 for unit cell geometry and associated meshing). As seen, each unit cell has a period of 400 nm and contains 24 gold pixels, each of which measures 50 nm × 50 nm × 20 nm. As before, we employ a magnetic current excitation with a sinusoidally modulated Gaussian amplitude defined in (18) for normal incidence.

Fig. 7. The geometry of a pixelized gold metasurface unit cell and its prismatic spatial discretization.

Download Full Size | PDF

Figure 8(a) shows the transient response of the metasurface when subject to an x-polarized normally incident illumination. Through Fourier transformation, reflection and transmission performance of the metasurface can be extracted as shown in Fig. 8(b). In the same way, the results from a y-polarized illumination are demonstrated in Fig. 9. In order to validate the accuracy and efficiency of the proposed method, the results are compared with the commercial CST software package. The results from the CST FEM solver with imported dispersion data are also plotted in Figs. 8(b) and 9(b). These results demonstrate good agreement between the CST FEM solver and the proposed algorithm. The prismatic DGTD + GDM method requires 312 seconds for the computation while CST takes 662 seconds. The efficiency improvement of the proposed method comes from the prismatic elements and the ability of the DGTD formulation to readily exploit parallel computing architectures.

Fig. 8. (a) Transient response of pixelized metasurface under x-polarized normally incident illumination. (b) Transmission and reflection performance in wavelengths.

Download Full Size | PDF

Fig. 9. (a) Transient response of pixelized metasurface under y-polarized normally incident illumination. (b) Transmission and reflection performance in wavelengths.

Download Full Size | PDF

The simulation toolkit has been made available online [41]. The user of the toolkit can follow the instructions, customize the pixelized gold pattern of the metasurface, and investigate the transient response as well as the transmission and reflection properties of the metasurface.

4. Conclusions

A prism DGTD algorithm together with the GDM model was proposed for simulating electromagnetic fields of planar photonic devices (e.g., FSS and metasurfaces). The ADE method and a high-order Runge-Kutta scheme were introduced as an effective methodology for integrating the GDM model into prism-based DGTD. The accuracy of the proposed algorithm was analyzed by comparing results with analytical the reflection and transmission properties of a thin gold film. The proposed algorithm was then used to simulate several numerical examples. The extracted S-parameters agree well with the results produced by commercial software packages. Moreover, the prism-based DGTG formulation was found to be more efficient than the conventional tetrahedral-based DGTD utilized in commercial solvers due to the prismatic spatial discretization with reduced unknowns for planar structures. This work introduces a new simulation algorithm for the efficient numerical analysis of planar photonic devices. A simulation toolkit with the proposed algorithm has been released online to public, enabling the users to efficiently analyze pixelized metasurfaces with customized patterns.

Funding

Penn State MRSEC; National Science Foundation (DMR-1420620); Defense Advanced Research Projects Agency (HR00111720032); Engineering and Physical Sciences Research Council (EP/R013918/1); National Natural Science Foundation of China (61901087, 61901132).

Acknowledgments

The authors would like to thank Dr. Huaguang Bao in Nanjing University of Science and Technology for discussions during the early stages of this work.

Disclosures

The authors declare no conflicts of interest.

References

1. A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, “Planar photonics with metasurfaces,” Science 339(6125), 1232009 (2013). [CrossRef]

2. P. Berto, L. Philippet, J. Osmond, C. F. Liu, A. Afridi, M. Montagut Marques, B. Molero Agudo, G. Tessier, and R. Quidant, “Tunable and free-form planar optics,” Nat. Photonics 13(9), 649–656 (2019). [CrossRef]

3. L. Bi, R. Khorasaninejad, B. Wessels, B. Stadler, and V. J. Sorger, “Planar oxide photonic materials and devices: feature issue introduction,” Opt. Mater. Express 8(10), 3250 (2018). [CrossRef]

4. J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Texts in Applied Mathematics No. 54 (Springer, 2008).

5. J. Alvarez, L. D. Angulo, M. F. Pantoja, A. R. Bretones, and S. G. Garcia, “Source and boundary implementation in vector and scalar DGTD,” IEEE Trans. Antennas Propag. 58(6), 1997–2003 (2010). [CrossRef]

6. J. Alvarez, L. D. Angulo, M. R. Cabello, A. R. Bretones, and S. G. Garcia, “An analysis of the leap-frog discontinuous Galerkin method for Maxwell’s equations,” IEEE Trans. Microwave Theory Tech. 62(2), 197–207 (2014). [CrossRef]

7. H. M. Lee, S. Gao, E. X. Liu, G. S. Samudra, and E. P. Li, “Two-dimensional discontinuous Galerkin time-domain method for modeling of arbitrarily shaped power-ground planes,” IEEE Trans. Electromagn. Compat. 57(6), 1744–1747 (2015). [CrossRef]

8. J. Alvarez, L. Angulo, A. Rubio Bretones, and S. G. Garcia, “A spurious-free discontinuous Galerkin time-domain method for the accurate modeling of microwave filters,” IEEE Trans. Microwave Theory Tech. 60(8), 2359–2369 (2012). [CrossRef]

9. R. Holland, V. P. Cable, and L. C. Wilson, “Finite-volume time-domain (FVTD) techniques for EM scattering,” IEEE Trans. Electromagn. Compat. 33(4), 281–294 (1991). [CrossRef]

10. J.-M. Jin, The Finite Element Method in Electromagnetics, 3rd edition (John Wiley & Sons Inc, 2014).

11. S. G. Garcia, R. G. Rubio, A. R. Bretones, and R. G. Martin, “Extension of the ADI-FDTD method to Debye media,” IEEE Trans. Antennas Propag. 51(11), 3183–3186 (2003). [CrossRef]

12. P. Li, L. J. Jiang, and H. Bağci, “Transient analysis of dispersive power-ground plate pairs with arbitrarily shaped antipads by the DGTD method with wave port excitation,” IEEE Trans. Electromagn. Compat. 59(1), 172–183 (2017). [CrossRef]

13. S. D. Gedney, J. C. Young, T. C. Kramer, and J. A. Roden, “A discontinuous Galerkin finite element time-domain method modeling of dispersive media,” IEEE Trans. Antennas Propag. 60(4), 1969–1977 (2012). [CrossRef]

14. W. K. Anderson, L. Wang, S. Kapadia, C. Tanis, and B. Hilbert, “Petrov–Galerkin and discontinuous-Galerkin methods for time-domain and frequency-domain electromagnetic simulations,” J. Comput. Phys. 230(23), 8360–8385 (2011). [CrossRef]

15. T. Lu, P. Zhang, and W. Cai, “Discontinuous Galerkin methods for dispersive and lossy Maxwell's equations and PML boundary conditions,” J. Comput. Phys. 200(2), 549–580 (2004). [CrossRef]

16. N. Schmitt, C. Scheid, S. Lanteri, A. Moreau, and J. Viquerat, “A DGTD method for the numerical modeling of the interaction of light with nanometer scale metallic structures taking into account non-local dispersion effects,” J. Comput. Phys. 316, 396–415 (2016). [CrossRef]

17. R. Léger, J. Viquerat, C. Durochat, C. Scheid, and S. Lanteri, “A parallel non-conforming multi-element DGTD method for the simulation of electromagnetic wave interaction with metallic nanoparticles,” J. Computat. Appl. Mathematics 270, 330–342 (2014). [CrossRef]

18. S.-C. Kong, J. J. Simpson, and V. Backman, “ADE-FDTD scattered-field formulation for dispersive materials,” IEEE Microw. Wireless Compon. Lett. 18(1), 4–6 (2008). [CrossRef]

19. A. Bruner, D. Eger, M. B. Oron, P. Blau, M. Katz, and S. Ruschin, “Temperature-dependent Sellmeier equation for the refractive index of stoichiometric lithium tantalate,” Opt. Lett. 28(3), 194–196 (2003). [CrossRef]

20. L. J. Prokopeva, J. D. Borneman, and A. V. Kildishev, “Optical dispersion models for time-domain modeling of metal-dielectric nanostructures,” IEEE Trans. Magn. 47(5), 1150–1153 (2011). [CrossRef]

21. Q. Ren, H. Bao, S. D. Campbell, L. J. Prokopeva, A. V. Kildishev, and D. H. Werner, “Continuous-discontinuous Galerkin time domain (CDGTD) method with generalized dispersive material (GDM) model for computational photonics,” Opt. Express 26(22), 29005 (2018). [CrossRef]

22. H. Bao, L. Kang, S. D. Campbell, and D. H. Werner, “Discontinuous Galerkin time domain analysis of electromagnetic scattering from dispersive periodic nanostructures at oblique incidence,” Opt. Express 27(9), 13116 (2019). [CrossRef]

23. M. L. N. Chen, L. J. Jiang, and W. E. I. Sha, “Artificial perfect electric conductor-perfect magnetic conductor anisotropic metasurface for generating orbital angular momentum of microwave with nearly perfect conversion efficiency,” J. Appl. Phys. 119(6), 064506 (2016). [CrossRef]

24. X. Yang, W. Zang, M. Ma, P. Zhan, Z. Chen, and Z. Wang, “Continuous phase control of second harmonic generation from metasurfaces composed of complementary split ring resonators,” Opt. Express 25(23), 28363 (2017). [CrossRef]

25. P.-C. Li, Y. Zhao, A. Alù, and E. T. Yu, “Experimental realization and modeling of a subwavelength frequency-selective plasmonic metasurface,” Appl. Phys. Lett. 99(22), 221106 (2011). [CrossRef]

26. W. Mai, D. Zhu, Z. Gong, X. Lin, Y. Chen, J. Hu, and D. H. Werner, “Broadband transparent chiral mirrors: Design methodology and bandwidth analysis,” AIP Adv. 9(4), 045305 (2019). [CrossRef]

27. I. Hussain, H. Li, and Q. Cao, “Multiscale structure simulation using adaptive mesh in DGTD method,” IEEE J. Multiscale Multiphys. Comput. Tech. 2, 115–123 (2017). [CrossRef]

28. W. Mai, J. Hu, P. Li, and H. Zhao, “An efficient and stable 2-D/3-D hybrid discontinuous Galerkin time-domain analysis with adaptive criterion for arbitrarily shaped antipads in dispersive parallel-plate pair,” IEEE Trans. Microwave Theory Tech. 65(10), 3671–3681 (2017). [CrossRef]

29. W. Mai, P. Li, C. G. Li, M. Jiang, W. Hao, L. Jiang, and J. Hu, “A straightforward updating criterion for 2-D/3-D hybrid discontinuous Galerkin time-domain method controlling comparative error,” IEEE Trans. Microwave Theory Tech. 66(4), 1713–1722 (2018). [CrossRef]

30. W. Mai, P. Li, H. Bao, X. Li, L. Jiang, J. Hu, and D. H. Werner, “Prism-based DGTD with a simplified periodic boundary condition to analyze FSS with D2n symmetry in a rectangular array under normal incidence,” Antennas Wirel. Propag. Lett. 18(4), 771–775 (2019). [CrossRef]

31. L. J. Prokopeva and A. V. Kildishev, “Efficient time-domain model of the graphene dielectric function,” in Proceedings of SPIE - The International Society for Optical Engineering (2013).

32. R. D. Graglia, D. R. Wilton, A. F. Peterson, and I. L. Gheorma, “Higher order interpolatory vector bases on prism elements,” IEEE Trans. Antennas Propag. 46(3), 442–450 (1998). [CrossRef]

33. P. Dular, J. Y. Hody, A. Nicolet, A. Genon, and W. Legros, “Mixed finite elements associated with a collection of tetrahedra, hexahedra and prisms,” IEEE Trans. Magn. 30(5), 2980–2983 (1994). [CrossRef]

34. L. D. Angulo, J. Alvarez, M. F. Pantoja, S. G. Garcia, and A. R. Bretones, “Discontinuous Galerkin time domain methods in computational electrodynamics: state of the art,” in Forum for Electromagnetic Research Methods and Application Technologies (FERMAT) (2014).

35. P. G. Etchegoin, E. C. Le Ru, and M. Meyer, “An analytic model for the optical properties of gold,” J. Chem. Phys. 125(16), 164705 (2006). [CrossRef]

36. P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B 6(12), 4370–4379 (1972). [CrossRef]

37. E. X. Jin and X. Xu, “Plasmonic effects in near-field optical transmission enhancement through a single bowtie-shaped aperture,” Appl. Phys. B 84(1-2), 3–9 (2006). [CrossRef]

38. D. H. Werner, Broadband Metamaterials in Electromagnetics (Pan Stanford Publishing, 2017).

39. D. Z. Zhu, M. D. Gregory, P. L. Werner, and D. H. Werner, “Fabrication and characterization of multiband polarization independent 3-D-printed frequency selective structures with ultrawide fields of view,” IEEE Trans. Antennas Propag. 66(11), 6096–6105 (2018). [CrossRef]

40. R. L. Haupt and D. H. Werner, Genetic Algorithms in Electromagnetics (Wiley-IEEE Press, 2007).

41. W. Mai and D. Werner, “Prism DGTD simulation toolkit with GDM for pixelized metasurfaces,” https://github.com/maiwending/pDGTD_GDM (2020).

Parameters	Drude Term	1^st Critical Point Term	2^nd Critical Point Term
$a_{0}$	$1.8103 \times 10^{32}$	$2.6976 \times 10^{31}$	$1.0587 \times 10^{32}$
$a_{1}$	0	$2.1508 \times 10^{15}$	$1.0239 \times 10^{16}$
$b_{0}$	0	$1.8534 \times 10^{31}$	$3.7523 \times 10^{31}$
$b_{1}$	$1.108 \times 10^{14}$	$1.5697 \times 10^{15}$	$3.9656 \times 10^{15}$

Method	Number of Elements^b	Number of Unknowns	CPU Time(s)	Memory (MB)
Tetra CST FEM with Exp Data	$234, 147$	$4, 400, 326$	$2658 s$	4,432
Tetra DGTD + GDM	285,831	$13, 719, 888$	$870 s$	15,456
Prism DGTD + GDM	47,340	$3, 408, 480$	$230 s$	5,126

Parameters	Drude Term	1^st Critical Point Term	2^nd Critical Point Term
$a_{0}$	$1.8103 \times 10^{32}$	$2.6976 \times 10^{31}$	$1.0587 \times 10^{32}$
$a_{1}$	0	$2.1508 \times 10^{15}$	$1.0239 \times 10^{16}$
$b_{0}$	0	$1.8534 \times 10^{31}$	$3.7523 \times 10^{31}$
$b_{1}$	$1.108 \times 10^{14}$	$1.5697 \times 10^{15}$	$3.9656 \times 10^{15}$

Method	Number of Elements^b	Number of Unknowns	CPU Time(s)	Memory (MB)
Tetra CST FEM with Exp Data	$234, 147$	$4, 400, 326$	$2658 s$	4,432
Tetra DGTD + GDM	285,831	$13, 719, 888$	$870 s$	15,456
Prism DGTD + GDM	47,340	$3, 408, 480$	$230 s$	5,126

Prismatic discontinuous Galerkin time domain method with an integrated generalized dispersion model for efficient optical metasurface analysis

Abstract

1. Introduction

2. Methods and formulation

2.1 Governing equations for the generalized dispersive model

2.2 Spatial discretization of the prism DGTD

2.3 Runge-Kutta time stepping scheme

2.4 Extraction of transmission and reflection

3. Numerical examples

3.1 Optimal settings of the generalized dispersive model

3.2 Reflection and transmission of a thin gold film

3.3 Illumination on a thin gold hole array

3.4 Pixelized optical metasurface

4. Conclusions

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (9)

Tables (2)

Equations (37)

Optical Materials Express