Discontinuous Galerkin time domain analysis of electromagnetic scattering from dispersive periodic nanostructures at oblique incidence

Huaguang Bao; Lei Kang; Sawyer D. Campbell; Douglas H. Werner

doi:10.1364/OE.27.013116

1. Introduction

The discontinuous Galerkin time domain (DGTD) method provides advanced features such as non-conformal meshing [1], element-level domain decomposition [2], and a high degree of parallelization [3] which has helped it attract considerable attention in recent years for its potential in flexible geometric modeling and fast electromagnetic simulations [4–6]. In the terahertz, infrared, and optical frequency regimes metals cannot be simply treated as perfect electric conductors as they are at radio frequencies. Therefore, an accurate and efficient dispersive model for metals and other dispersive materials is a subject of great importance in computational photonics. For noble metals and plasmonic materials, the Drude-critical point (DCP) model has shown its own superiority over Drude and Drude-Lorentz models with physics-based dispersion approximations [7,8]. For dispersive dielectrics and semiconductor alloys, the Debye, Tauc-Lorentz, generalized oscillator, and generalized critical points models along with classical empirical models have been introduced to represent the material dispersion [9,10]. A generalized dispersive material (GDM) model based on [1/2] Padé approximants which covers Debye, Drude, Lorentz, Sellmeier, and critical points models with an arbitrary number of poles has been proposed to achieve a robust universal model for various dispersive materials [11–13]. The GDM model has been successfully incorporated into the DGTD method and its implementation with advanced basis functions enables more accurate simulation with extremely coarse meshes compared to conventional finite element time domain (FETD) and finite difference time domain (FDTD) methods [14].

Periodic nanostructures remain a high priority in many photonic material, metasurface, and metamaterial applications [15–17]. To reduce the computational complexity and enable fast analysis of these periodic nanostructures, their periodic nature should be fully exploited by reducing the infinite structure to a single unit cell with appropriate boundary conditions. In the general case of oblique incidence, the interaction between incident plane waves and periodic structures can be modeled by Floquet periodicity. For time domain solvers, the time shift between the tangential field components on the parallel periodic surfaces makes the implementation of periodic boundary conditions (PBCs) not straightforward. The constant wavenumber approach [18–20] and the transformed field variable technique [21–23] have been proposed as two principle means to embed this time shift. As a spectral technique, the former is directly derived from the standard Maxwell’s equations, which preserves the formulations for dispersive structures [19] and absorbing boundary conditions [20]. However, it does not provide a solution for time-dependent fields, which makes it inapplicable in some multi-physics and nonlinear applications [24,25]. In contrast, the latter generates transformed governing equations different from the standard Maxwell’s equations, thus the relevant formulations need to be rederived. Nevertheless, the time-dependent fields can be obtained by a reverse transformation.

In this paper, the DGTD method with a GDM model has been extended to the electromagnetic analysis of periodic nanostructures under obliquely incident illumination. The transformed field variable technique is introduced to implement the PBCs applicable for both normal and oblique incidence cases. Based on the transformed Maxwell’s equations, the formulations of the proposed DGTD method with a GDM model are derived using an auxiliary differential equation (ADE) technique. An s–stage Runge-Kutta (RK) time integration scheme is proposed to retain a high order accuracy of the DGTD method. The proposed method is validated by analyzing several periodic nanostructures with dispersive elements, including reflection and transmission of a thin film, surface plasmon excitation at the interfaces of a metallic hole array structure, and absorption of a dual-band infrared absorber. This work reveals a possible potential of DGTD-based techniques to multiscale and multi-physics design as well as optimization of periodic nanostructures.

2. Methods and formulations

2.1. Transformed Maxwell’s equations and PBCs

When a doubly periodic structure is illuminated by a time-harmonic plane wave at an oblique incident angle (θ_i, φ_i) as illustrated in Fig. 1, the fields on the left/bottom periodic cell boundary can be obtained from the fields on the right/top boundary by multiplying it with a corresponding Floquet phase factor [22]

\vec{E} (x, y, z, ω) = {\begin{matrix} \vec{E} (x + x_{p}, y, z, ω) e^{j k_{0} x_{p} \sin θ_{i} \cos φ_{i}} \\ \vec{E} (x, y + y_{p}, z, ω) e^{j k_{0} y_{p} \sin θ_{i} \cos φ_{i}} \end{matrix}

\vec{H} (x, y, z, ω) = {\begin{matrix} \vec{H} (x + x_{p}, y, z, ω) e^{j k_{0} x_{p} \sin θ_{i} \cos φ_{i}} \\ \vec{H} (x, y + y_{p}, z, ω) e^{j k_{0} y_{p} \sin θ_{i} \cos φ_{i}} \end{matrix}

where x_p and y_p are the period lengths in the x- and y-directions, respectively, and k₀ is the wavenumber in free space. In the case of oblique incidence, these phase shifts in Eqs. (1) and (2) are translated from spatial shifts in the frequency domain to temporal shifts in the time domain. Here, the following auxiliary variables are introduced to embed these shifts

Fig. 1 Illustration of a doubly periodic structure under oblique illumination.

Download Full Size | PDF

{\vec{P}}_{e} (x, y, z, ω) = \vec{E} (x, y, z, ω) e^{j (k_{0} x \sin θ_{i} \cos φ_{i} + k_{0} y \sin θ_{i} \sin φ_{i})}

{\vec{P}}_{h} (x, y, z, ω) = \vec{H} (x, y, z, ω) e^{j (k_{0} x \sin θ_{i} \cos φ_{i} + k_{0} y \sin θ_{i} \sin φ_{i})}

With the transformed field variables, the periodic boundary conditions of Eqs. (1) and (2) can be reformulated as

{\vec{P}}_{e} (x, y, z, ω) = {\begin{matrix} {\vec{P}}_{e} (x + x_{p}, y, z, ω) \\ {\vec{P}}_{e} (x, y + y_{p}, z, ω) \end{matrix}

{\vec{P}}_{h} (x, y, z, ω) = {\begin{matrix} {\vec{P}}_{h} (x + x_{p}, y, z, ω) \\ {\vec{P}}_{h} (x, y + y_{p}, z, ω) \end{matrix} .

Substituting Eqs. (3) and (4) into the first order Maxwell’s equations in the frequency domain, the following transformed Maxwell’s equations are obtained:

j ω ε_{0} ε_{r} {\vec{P}}_{e} = \nabla \times {\vec{P}}_{h} - j ω \frac{1}{c} {\vec{k}}_{t} \times {\vec{P}}_{h}

j ω μ_{0} μ_{r} {\vec{P}}_{h} = - \nabla \times {\vec{P}}_{e} + j ω \frac{1}{c} {\vec{k}}_{t} \times {\vec{P}}_{e}

where

{\vec{k}}_{t} = (\sin θ_{i} \cos φ_{i}, \sin θ_{i} \sin φ_{i}, 0)

defines the transverse component of the unit incident vector wavenumber and c is speed of light in free space.

2.2. Generalized dispersive material model

The GDM model [11–14] represents classical dispersion formulas such as the Debye, Drude, Lorentz, Sellmeier, and critical points models as particular cases of a [1/2] Padé approximant. The relative permittivity of dispersive media described with the GDM model can be expressed as a sum of a finite constant permittivity $ε_{\infty}$ and multiple $χ_{m}$ terms

ε_{r} (ω) = ε_{\infty} + \sum_{m} χ_{m} (ω) = ε_{\infty} + \sum_{m} \frac{a_{0, m} + j ω a_{1, m}}{b_{0, m} + j ω b_{1, m} - ω^{2}} .

The values of GDM parameters a₀, a₁, b₀, and b₁ for different dispersion models are shown in Table 1, while the parameters and their physical meanings in Table 1 can be found in [11,14].

Table 1. GDM Parameters for Different Dispersion Terms

View Table | View all tables in this article

By substituting Eq. (9) into Eq. (7) and introducing the following auxiliary differential equation

{\vec{Q}}_{m} = ε_{0} \frac{a_{0, m} + j ω a_{1, m}}{b_{0, m} + j ω b_{1, m} - ω^{2}} {\vec{P}}_{e},

the transformed Maxwell’s equations with a GDM model can be expressed as

j ω ε_{0} ε_{\infty} {\vec{P}}_{e} + \sum_{m} j ω {\vec{Q}}_{m} = \nabla \times {\vec{P}}_{h} - j ω \frac{1}{c} {\vec{k}}_{t} \times {\vec{P}}_{h}

Transforming Eqs. (5), (6), (8), (10), and (11) into the time domain, the transformed Maxwell’s equations with a GDM model and PBCs in the time domain can be obtained as

ε_{0} ε_{\infty} \frac{\partial {\vec{P}}_{e}}{\partial t} + \sum_{m} \frac{\partial {\vec{Q}}_{m}}{\partial t} = \nabla \times {\vec{P}}_{h} - \frac{1}{c} {\vec{k}}_{t} \times \frac{\partial {\vec{P}}_{h}}{\partial t}

μ_{0} μ_{r} \frac{\partial {\vec{P}}_{h}}{\partial t} = - \nabla \times {\vec{P}}_{e} + \frac{1}{c} {\vec{k}}_{t} \times \frac{\partial {\vec{P}}_{e}}{\partial t}

\frac{\partial^{2} {\vec{Q}}_{m}}{\partial t^{2}} + b_{1, m} \frac{\partial {\vec{Q}}_{m}}{\partial t} + b_{0, m} {\vec{Q}}_{m} = ε_{0} a_{0, m} {\vec{P}}_{e} + ε_{0} a_{1, m} \frac{\partial {\vec{P}}_{e}}{\partial t}

{\vec{P}}_{e} (x, y, z, t) = {\begin{matrix} {\vec{P}}_{e} (x + x_{p}, y, z, t) \\ {\vec{P}}_{e} (x, y + y_{p}, z, t) \end{matrix}

{\vec{P}}_{h} (x, y, z, t) = {\begin{matrix} {\vec{P}}_{h} (x + x_{p}, y, z, t) \\ {\vec{P}}_{h} (x, y + y_{p}, z, t) \end{matrix}

It should be stated that in this work, the computational domain in the non-periodic directions is truncated by an anisotropic perfectly matched layer (PML) absorbing material [26,27].

2.3. Discontinuous Galerkin spatial discretization

As in standard DGTD analysis, the entire computational domain is first discretized with smaller tetrahedral or hexahedral elements. Next, the variables in the transformed Maxwell’s equations with the GDM model can be expanded in terms of hierarchical vector basis functions $\vec{Φ}$ [28,29]:

{\vec{P}}_{e} = \sum p_{e} \vec{Φ}, {\vec{P}}_{h} = \sum p_{h} \vec{Φ}, {\vec{Q}}_{m} = \sum q_{m} \vec{Φ}

where p_e, p_h, and q_m are the unknown expansion coefficients. Using a weighted residual method referred to as Galerkin testing, the transformed Maxwell’s equations are tested by

\vec{Φ}

and reformulated as

\begin{array}{l} ∭ \vec{Φ} \cdot ε_{0} ε_{\infty} \frac{\partial {\vec{P}}_{e}}{\partial t} d V + \sum_{m} ∭ \vec{Φ} \cdot \frac{\partial {\vec{Q}}_{m}}{\partial t} d V \\ = ∭ \nabla \times \vec{Φ} \cdot {\vec{P}}_{h} d V - ∭ \vec{Φ} \cdot (\frac{1}{c} {\vec{k}}_{t} \times \frac{\partial {\vec{P}}_{h}}{\partial t}) d V + \underset{s}{∯} \vec{Φ} \cdot (\vec{n} \times {\vec{P}}_{h}) d S \end{array}

∭ \vec{Φ} \cdot μ_{0} μ_{r} \frac{\partial {\vec{P}}_{h}}{\partial t} d V = - ∭ \nabla \times \vec{Φ} \cdot {\vec{P}}_{e} d V + ∭ \vec{Φ} \cdot (\frac{1}{c} {\vec{k}}_{t} \times \frac{\partial {\vec{P}}_{e}}{\partial t}) d V - \underset{s}{∯} \vec{Φ} \cdot (\vec{n} \times {\vec{P}}_{e}) d S

\begin{array}{l} ∭ \vec{Φ} \cdot \frac{\partial^{2} {\vec{Q}}_{m}}{\partial t^{2}} d V + b_{1, m} ∭ \vec{Φ} \cdot \frac{\partial {\vec{Q}}_{m}}{\partial t} d V + b_{0, m} ∭ \vec{Φ} \cdot {\vec{Q}}_{m} d V \\ = ε_{0} a_{0, m} ∭ \vec{Φ} \cdot {\vec{P}}_{e} d V + ε_{0} a_{1, m} ∭ \vec{Φ} \cdot \frac{\partial {\vec{P}}_{e}}{\partial t} d V . \end{array}

The tangential continuity of the electrical and magnetic fields at the interfaces between two adjacent elements is imposed by modifying the surface integral terms in Eqs. (18) and (19) with an appropriate numerical flux. Here, the upwind flux [23] is employed by taking the following forms with the transformed field variables:

{\begin{matrix} \vec{n} \times {\vec{P}}_{e} = \vec{n} \times \frac{(Y {\vec{P}}_{e} - \vec{n} \times {\vec{P}}_{h}) + (Y^{+} {\vec{P}}_{e}^{+} - \vec{n} \times {\vec{P}}_{h}^{+})}{Y + Y^{+}} \\ \vec{n} \times {\vec{P}}_{h} = \vec{n} \times \frac{(Z {\vec{P}}_{h} - \vec{n} \times {\vec{P}}_{e}) + (Z^{+} {\vec{P}}_{h}^{+} - \vec{n} \times {\vec{P}}_{e}^{+})}{Z + Z^{+}} \end{matrix}

where

Z = 1 / Y

is the intrinsic impedance and the superscript

" + "

implies the variables from an adjacent element. Particularly, for the fields on a specific periodic boundary, the superscript

" + "

signifies the variables from the corresponding element on the opposite parallel periodic surface.

By substituting Eqs. (17) and (21) into Eqs. (18)–(20), the semi-discrete DGTD matrix equations with the GDM model can be obtained

T_{e e}^{ε_{0} ε_{\infty}} \frac{\partial p_{e}}{\partial t} + M_{e h} \frac{\partial p_{h}}{\partial t} + T_{e e}^{1} \frac{\partial q_{m}}{\partial t} = S_{e e} p_{e} + S_{e h} p_{h} + F_{e e} p_{e}^{+} + F_{e h} p_{h}^{+}

T_{h h}^{μ_{0} μ_{r}} \frac{\partial p_{h}}{\partial t} - M_{h e} \frac{\partial p_{e}}{\partial t} = S_{h h} p_{h} + S_{h e} p_{e} + F_{h h} p_{h}^{+} + F_{h e} p_{e}^{+}

\frac{\partial^{2} q_{m}}{\partial t^{2}} + b_{1, m} \frac{\partial q_{m}}{\partial t} - ε_{0} a_{1, m} \frac{\partial p_{e}}{\partial t} = ε_{0} a_{0, m} p_{e} - b_{0, m} q_{m}

where p_e, p_h, and q_m are the unknown vectors and T, M, and S are sparse matrices for which the subscripts indicate the corresponding test and basis functions. The matrix elements are defined as

\begin{matrix} {[T_{r s}^{ξ}]}_{i j} = ξ ∭ {\vec{Φ}}_{r i} \cdot {\vec{Φ}}_{s j} d V \\ {[M_{r s}]}_{i j} = ∭ {\vec{Φ}}_{r i} \cdot (\frac{1}{c} {\vec{k}}_{t} \times {\vec{Φ}}_{s j}) d V \\ {[S_{e e}]}_{i j} = - ∯_{s} {\vec{Φ}}_{e i} \cdot \frac{\vec{n} \times \vec{n} \times {\vec{Φ}}_{e j}}{Z + Z^{+}} d S \\ {[S_{e h}]}_{i j} = ∭ \nabla \times {\vec{Φ}}_{e i} \cdot {\vec{Φ}}_{h j} d V + ∯_{s} {\vec{Φ}}_{e i} \cdot \frac{\vec{n} \times Z {\vec{Φ}}_{h j}}{Z + Z^{+}} d S \\ {[F_{e e}]}_{i j} = - ∯_{s} {\vec{Φ}}_{e i} \cdot \frac{\vec{n} \times \vec{n} \times {\vec{Φ}}_{e j}^{+}}{Z + Z^{+}} d S \\ {[F_{e h}]}_{i j} = ∯_{s} {\vec{Φ}}_{e i} \cdot \frac{\vec{n} \times Z^{+} {\vec{Φ}}_{h j}^{+}}{Z + Z^{+}} d S \\ {[S_{h h}]}_{i j} = ∯_{s} {\vec{Φ}}_{h i} \cdot \frac{\vec{n} \times \vec{n} \times {\vec{Φ}}_{h j}}{Y + Y^{+}} d S \\ {[S_{h e}]}_{i j} = - ∭ \nabla \times {\vec{Φ}}_{h i} \cdot {\vec{Φ}}_{e j} d V - ∯_{s} {\vec{Φ}}_{h i} \cdot \frac{\vec{n} \times Y {\vec{Φ}}_{h j}}{Y + Y^{+}} d S \\ {[F_{h h}]}_{i j} = ∯_{s} {\vec{Φ}}_{h i} \cdot \frac{\vec{n} \times \vec{n} \times {\vec{Φ}}_{h j}^{+}}{Y + Y^{+}} d S \\ {[F_{h e}]}_{i j} = - ∯_{s} {\vec{Φ}}_{h i} \cdot \frac{\vec{n} \times Y^{+} {\vec{Φ}}_{e j}^{+}}{Y + Y^{+}} d S \end{matrix}

where

r = e, h

,

s = e, h

, and

ξ = 1, ε_{0} ε_{\infty}, μ_{0} μ_{r}

.

2.4. Runge-Kutta time stepping scheme

Due to the second-order time derivative term in Eq. (24), the discontinuous Galerkin semi-discrete system of the transformed GDM Maxwell’s equations should be remodified to make an s–stage Runge-Kutta time integration scheme available. Here, the following auxiliary differential equation is introduced:

\frac{\partial q_{m}}{\partial t} = γ_{m} .

With the auxiliary differential equation, the semi-discrete DGTD matrix equations can be reformulated as

M \frac{\partial u}{\partial t} = L_{1} u + L_{1} u^{+}

where u = [ p_e, p_h, q_m, γ_m]^T is the unknown vector and the associated matrix quantities are given by

M = [\begin{matrix} T_{e e}^{ε_{0} ε_{\infty}} & M_{e h} & 0 & 0 \\ - M_{h e} & T_{h h}^{μ_{0} μ_{r}} & 0 & 0 \\ - ε_{0} a_{1, m} I & 0 & 0 & I \\ 0 & 0 & I & 0 \end{matrix}]

L_{1} = [\begin{matrix} S_{e e} & S_{e h} & 0 & - T_{e q}^{1} \\ S_{h e} & S_{h h} & 0 & 0 \\ ε_{0} a_{0, m} I & 0 & - b_{0, m} I & - b_{1, m} I \\ 0 & 0 & 0 & I \end{matrix}]

L_{2} = [\begin{matrix} F_{e e} & F_{e h} & 0 & 0 \\ F_{h e} & F_{h h} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}]

where I is the identity matrix. At this point using an s–stage RK time integration scheme, the fully discrete DGTD system is obtained as

{\begin{matrix} u (t + c_{i} Δ t) = u (t) + Δ t \sum_{j = 1}^{s} a_{i j} M^{- 1} (L_{1} u (t + c_{j} Δ t) + L_{1} u^{+} (t + c_{j} Δ t)), 1 \leq i \leq s \\ u (t + Δ t) = u (t) + Δ t \sum_{i = 1}^{s} b_{i} M^{- 1} (L_{1} u (t + c_{i} Δ t) + L_{1} u^{+} (t + c_{i} Δ t)) \end{matrix}

where

Δ t

is the time step increment and the fixed scalar coefficients

a_{i j}

,

b_{i}

, and

c_{i}

can be found in [30].

3. Numerical results

In this section, a Drude-critical points model composed of one Drude term and two critical points terms (D2CP) is employed to describe the dispersion of gold, which has demonstrated improved accuracy and better efficiency than the conventional Drude-Lorentz model. The detailed parameters of the D2CP model have been given in [14]. For all of the following examples, two PMLs are used to truncate the computational domain in the z-direction and PBCs are assigned in the lateral direction. In order to obtain a wide-band response, a modulated Gaussian pulse is employed to generate the incident field, which is specified as

{\vec{E}}^{i n c} (z, t) = e^{- {(t - \frac{\cos θ_{i} (z - z_{s})}{t_{w}} - t_{s})}^{2}} \sin [2 π f_{0} (t - \frac{\cos θ_{i} (z - z_{s})}{t_{w}} - t_{s})]

where f₀ is the carrier frequency, t_w controls the bandwidth of the pulse, while z_s and t_s are the spatial and time shifts, respectively. For all cases, the incident electric field is assumed to be parallel to the x-axis and a fixed horizontal azimuth φ_i = 90° is used. The order number of spatial basis functions is two and the fields are updated in a low-storage fourth order RK scheme. As an explicit time domain solver, the RK-DGTD method is conditionally stable and its time increment should obey a specific stability condition such as the CFL stability criterion in the FDTD method. Details regarding the stability analysis of the DGTD method can be found in [31,32]. Generally, the maximum value of the time increment depends on the geometrical sizes of the elements and the propagation speed of electromagnetic waves in the media. Here, the time step increment is determined by the following estimate

Δ t \leq \frac{d_{m i n} \sqrt{ε_{r, \infty}}}{5 c}

where d_min is the minimal value of the diameter of an inscribed sphere for each element. Particularly, this time step increment decreases with an increasing incident pitch angle from normal incidence to grazing incidence.

3.1. Reflection and transmission of a thin film

An analytical canonical case consisting of a 50 nm thick gold film illuminated by a plane wave at different incident angles is simulated to validate the proposed DGTD-GDM-PBCs method in analyzing both normally and obliquely illuminated periodic nanostructures. The simulation setup is illustrated in Fig. 2(a). The investigated wavelength range extends from 200 nm to 2000 nm with a modulated Gaussian pulse whose detailed parameters have been provided in [14]. For normally incident plane wave illumination, the periodic boundaries can be defined as pairs of perfect electric conducting (PEC) and perfect magnetic conducting (PMC) boundary conditions [14] to serve as a comparison to the general PBC method proposed here. The transient electric fields at the probes are shown in Fig. 2(b) and the results agree well with those given in [14]. Further, the reflection and transmission coefficients under plane wave illumination with different pitch angles θ_i are calculated using a fast Fourier transform and compared to the analytical solutions. Figure 2(c)–2(e) demonstrates good agreement between the results from the proposed method and the analytical solutions, validating the proposed method under both normal and oblique incidence conditions.

Fig. 2 Gold film with D2CP model under plane wave illumination. (a) The simulation setup. (b) The transient reflected and transmitted E field at probes 1 and 2. (c,d,e) The reflection and transmission coefficients under plane wave illumination with different pitch angles (c) θ_i = 0°, (d) θ_i = 30°, and (e) θ_i = 60°.

Download Full Size | PDF

3.2. Surface plasmon excitation at the interface of a metal hole array

A nanohole array structure based on the 50-nanometer-thick gold film simulated above is generated by perforating 150 nm radius holes spaced with a 600 nm lattice period and evaluated using the proposed DGTD-GDM-PBCs method. As the unit cell schematic in Fig. 3(a) shows, the structure is illuminated by a normally incident plane wave. The reflection and transmission coefficients in the wavelength range of [600 nm, 2000 nm] are calculated by using both the proposed method and the CST software package [33] in which a frequency domain solver (finite element method, FEM) was utilized. As shown in Fig. 3(b), good agreement is demonstrated, confirming the effectiveness of the proposed DGTD-GDM-PBCs method for dispersive periodic nanostructure analysis. Upon further inspection, a remarkable transmission enhancement compared to the transmission spectrum of a bare film (Fig. 2(c)) is observed at wavelengths around 780 nm. As presented in [34], this type of extraordinary transmission behavior can be primarily attributed to the excitation of a surface plasmon at the metallic hole array structure interface. In brief, 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, resulting in unique features in the transmission and reflection spectra. To better understand the observed response of the gold hole array structure, a cut of the electric field distribution (z-component) at the resonance wavelength (783.5 nm) in the y = –300 nm plane is shown in Fig. 3(c), where strong surface plasmon excitation is observed. Moreover, a comparison of computational efficiency between the proposed method and the commercial software package (CST-FEM) is made here. Firstly, it should be noted that a transient analysis of periodic structures under obliquely incident illumination is not available in the current version of CST or COMSOL commercial software. Our previous work [14] has shown the superiority of the DGTD method with a GDM model, in which a convergence behavior beyond the conventional FETD and FDTD methods was exhibited. Therefore, a comparison of computational efficiency between the proposed method and the CST software package with a frequency domain FEM solver is executed on a DELL computer with an Intel(R) Xeon(R) E5-2680 CPU at 2.40 GHz and with 256 GB of RAM. As a time domain solver, the proposed method is able to obtain a wide-band response with a single simulation. The parameters of the incident pulse are f₀ = 350 THz, z_s = –375 nm, t_w = 2.0/(πf₀), and t_s = 3t_w for the frequency range of [100 THz, 700 THz]. The total number of time steps is chosen to be 2000 with a step length of 0.023 fs. In the CST software package, an adaptive frequency sweep tool is employed and a total of 33 frequency points are calculated to obtain the wide-band response with the same frequency range. The computing resources for the proposed method and the CST software package are summarized in Table 2.

Fig. 3 Surface plasmon excitation enabled extraordinary transmission of a metallic hole array structure. (a) Unit cell geometry. (b) The reflection and transmission coefficients at normal incidence. (c) Electric field distribution (E_z) in a cut plane at y = –300 nm.

Download Full Size | PDF

Table 2. Comparison on Computing Resources

View Table | View all tables in this article

Particularly, non-conformal curved hexahedrons are employed to discretize the entire computation domain in the proposed DGTD-GDM-PBCs method while conformal tetrahedrons are utilized in the CST software package. As observed in Table 2, the proposed method requires less computational time, which is reduced by about a factor of 2.7 when compared to the CST software package with the frequency domain FEM solver. The higher memory consumption of the proposed method is attributed to the non-conformal mesh and more matrix elements with the discontinuous Galerkin scheme. Furthermore, it should be noted that as a time domain solver, the proposed method is applicable to not only efficient multi-scale analysis, but also in multi-physics and nonlinear modeling.

3.3. Absorption of a dual-band infrared absorber

The performance of a dual-band infrared absorber [35,36] under plane wave illumination with a wide pitch angle range of [0°, 50°] is evaluated with the proposed method to further demonstrate its efficacy. The unit cell geometry of the absorber is illustrated in Fig. 4(a). The structure of the absorber consists of a 20-nanometer-thick gold elliptical nanodisk array, an optically thick (150 nm) gold substrate and a 30-nanometer-thick dielectric spacer (ε_r = 2.5). The reflectance spectra of the absorber under plane wave illumination with different pitch angles θ_i are calculated using the proposed DGTD method. In contrast, simulations of this configuration are executed with the CST software package using a frequency domain FEM solver. As shown in Fig. 4(b)–4(d), good agreement is demonstrated at both normal and oblique incidence. To demonstrate the performance of the absorber, the absorbance as a function of both wavelength and incident pitch angle is plotted in Fig. 4(e), which shows that the center frequencies of the two resonance peaks are maintained and the associated absorbance is higher than 83% for all the incident pitch angles up to 50 degrees.

Fig. 4 A dual-band infrared absorber. (a) Geometry of a unit cell. (b,c,d) Reflectance spectra under plane wave illuminations with different pitch angles (b) θ_i = 0°, (c) θ_i = 15°, and (d) θ_i = 30°. (e) Absorptivity as a function of both wavelength and pitch angle of incidence.

Download Full Size | PDF

4. Conclusions

In this paper, an efficient DGTD method with a generalized dispersive material model has been proposed for periodic dispersive structures analysis. The periodic boundary conditions were introduced using a transformed field variable technique. Based on the transformed Maxwell’s equations and with an ADE approach, the corresponding formulation of the DGTD method with a GDM model was derived. An s–stage Runge-Kutta time integration scheme has been utilized to update the semi-discrete equations by introducing an auxiliary differential equation while retaining the DGTD feature of high order accuracy. In addition, the proposed method has been successfully applied to obtain the wideband response of periodic nanostructures with dispersive elements at both normal and oblique incidence. The proposed framework is expected to provide an effective simulation tools for designing and optimizing novel nanostructured photonic elements.

Funding

DARPA/DSO Extreme Optics and Imaging (EXTREME) Program (HR00111720032).

References

1. S. Yan, C.-P. Lin, R. R. Arslanbekov, V. I. Kolobov, and J.-M. Jin, “A discontinuous Galerkin time-domain method with dynamically adaptive Cartesian mesh for computational electromagnetics,” IEEE. Trans. Antennas Propagat. 65(6), 3122–3133 (2017). [CrossRef]

2. H. Xu, D. Ding, J. Bi, and R. S. Chen, “A 3-D continuous-discontinuous Galerkin finite-element time-domain method for Maxwell’s equations,” IEEE Antennas Wirel. Propag. Lett. 16, 908–911 (2017). [CrossRef]

3. S. Dosopoulos, B. Zhao, and J. F. Lee, “Non-conformal and parallel discontinuous Galerkin time domain method for Maxwell’s equations: EM analysis of IC packages,” J. Comput. Phys. 238, 48–70 (2013). [CrossRef]

4. J. Alvarez, L. D. Angulo, A. R. Bretones, C. M. de Jong van Coevorden, and S. G. Garcia, “Efficient antenna modeling by DGTD: Leap-frog discontinuous Galerkin time-domain method,” IEEE Antennas Propag. Mag. 57(3), 95–106 (2015). [CrossRef]

5. J. Chen and Q. H. Liu, “Discontinuous Galerkin time domain methods for multiscale electromagnetic simulations: A review,” Proc. IEEE 101(2), 242–254 (2013). [CrossRef]

6. K. Stannigel, M. König, J. Niegemann, and K. Busch, “Discontinuous Galerkin time-domain computations of metallic nanostructures,” Opt. Express 17(17), 14934–14947 (2009). [CrossRef] [PubMed]

7. 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] [PubMed]

8. A. Vial and T. Laroche, “Comparison of gold and silver dispersion laws suitable for FDTD simulations,” Appl. Phys. B 93(1), 139–143 (2008). [CrossRef]

9. H. Chen and W. Z. Shen, “Perspectives in the characteristics and applications of Tauc-Lorentz dielectric function model,” Eur. Phys. J. B Cond. Matter Complex Syst. 43(4), 503–507 (2005). [CrossRef]

10. 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] [PubMed]

11. 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]

12. M. D. Thoreson, J. Fang, A. V. Kildishev, L. J. Prokopeva, P. Nyga, U. K. Chettiar, V. M. Shalaev, and V. P. Drachev, “Fabrication and realistic modeling of three-dimensional metal-dielectric composites,” J. Nanophotonics 5(1), 051513 (2011). [CrossRef]

13. N. K. Emani, D. Wang, T.-F. Chung, L. J. Prokopeva, A. V. Kildishev, V. M. Shalaev, Y. P. Chen, and A. Boltasseva, “Plasmon resonance in multilayer graphene nanoribbons,” Laser Photonics Rev. 9(6), 650–655 (2015). [CrossRef]

14. 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–29016 (2018). [CrossRef] [PubMed]

15. S. Y. Lin, J. G. Fleming, D. L. Hetherington, B. K. Smith, R. Biswas, K. M. Ho, M. M. Sigalas, W. Zubrzycki, S. R. Kurtz, and J. Bur, “A three-dimensional photonic crystal operating at infrared wavelengths,” Nature 394(6690), 251–253 (1998). [CrossRef]

16. H. T. Chen, A. J. Taylor, and N. Yu, “A review of metasurfaces: Physics and applications,” Rep. Prog. Phys. 79(7), 076401 (2016). [CrossRef] [PubMed]

17. S. Jeon, J.-U. Park, R. Cirelli, S. Yang, C. E. Heitzman, P. V. Braun, P. J. Kenis, and J. A. Rogers, “Fabricating complex three-dimensional nanostructures with high-resolution conformable phase masks,” Proc. Natl. Acad. Sci. U.S.A. 101(34), 12428–12433 (2004). [CrossRef] [PubMed]

18. A. Aminian and Y. Rahmat-Samii, “Spectral FDTD: A novel technique for the analysis of oblique incident plane wave on periodic structures,” IEEE Trans. Antenn. Propag. 54(6), 1818–1825 (2006). [CrossRef]

19. K. Elmahgoub, F. Yang, and A. Elsherbeni, Scattering Analysis of Periodic Structures Using Finite-Difference Time-Domain Method (San Rafael, CA, USA: Morgan & Claypool, 2012).

20. K. Sirenko, Y. Sirenko, and H. Bağcı, “Exact absorbing boundary conditions for periodic three-dimensional structures: Derivation and implementation in discontinuous Galerkin time-domain method,” IEEE J. Multiscale Multiphys. Comput. Techn. 3, 108–120 (2018). [CrossRef]

21. M. E. Veysoglu, R. T. Shin, and J. A. Kong, “A finite-difference time domain analysis of wave scattering from periodic structures: Oblique incidence case,” J. Electromagn. Waves Appl. 7(12), 1595–1607 (1993). [CrossRef]

22. L. E. R. Petersson and J. M. Jin, “A three-dimensional time-domain finite-element formulation for periodic structures,” IEEE Trans. Antenn. Propag. 54(1), 12–19 (2006). [CrossRef]

23. N. C. Miller, A. D. Baczewski, J. D. Albrecht, and B. Shanker, “A discontinuous Galerkin time domain framework for periodic structures subject to oblique excitation,” IEEE Trans. Antenn. Propag. 62(8), 4386–4391 (2014). [CrossRef]

24. K. Lopata and D. Neuhauser, “Multiscale Maxwell-Schrodinger modeling: A split field finite-difference time-domain approach to molecular nanopolaritonics,” J. Chem. Phys. 130(10), 104707 (2009). [CrossRef] [PubMed]

25. S. Yan, A. D. Greenwood, and J.-M. Jin, “Simulation of high-power microwave air breakdown modeled by a coupled Maxwell–Euler system with a non-Maxwellian EEDF,” IEEE. Trans. Antennas Propagat. 66(4), 1882–1893 (2018). [CrossRef]

26. S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE Trans. Antenn. Propag. 44(12), 1630–1639 (1996). [CrossRef]

27. H. Bao, L. Kang, S. D. Campbell, and D. H. Werner, “PML implementation in a non-conforming mixed-element DGTD method for periodic structure analysis,” (unpublished).

28. J. M. Jin, The Finite Element Method in Electromagnetics (Wiley, 2002).

29. M. M. Ilic and B. M. Notaros, “Higher order hierarchical curved hexahedral vector finite elements for electromagnetic modeling,” IEEE Trans. Antenn. Propag. 51(3), 1026–1033 (2003).

30. J. C. Butcher, Numerical Methods for Ordinary Differential Equations (Wiley, 2003).

31. L. Fezoui, S. Lanteri, S. Lohrengel, and S. Piperno, “Convergence and stability of a discontinuous Galerkin time-domain method for the 3-D heterogeneous Maxwell equations on unstructured meshes,” ESAIM: M2AN 39(6), 1149–1176 (2005). [CrossRef]

32. D. Sármány, M. A. Botchev, and J. J. W. van der Vegt, “Dispersion and dissipation error in high-order Runge–Kutta discontinuous Galerkin discretisations of the Maxwell equations,” J. Sci. Comput. 33(1), 47–74 (2007). [CrossRef]

33. C. S. T. Studio Suite, “CST Microwave Studio,” 2017. http://www.cst.com

34. C. Genet, M. P. van Exter, and J. P. Woerdman, “Fano-type interpretation of red shifts and red tails in hole array transmission spectra,” Opt. Express Commun. 225(4-6), 331–336 (2003). [CrossRef]

35. B. Zhang, Y. Zhao, Q. Hao, B. Kiraly, I.-C. Khoo, S. Chen, and T. J. Huang, “Polarization-independent dual-band infrared perfect absorber based on a metal-dielectric-metal elliptical nanodisk array,” Opt. Express 19(16), 15221–15228 (2011). [CrossRef] [PubMed]

36. Q. Ren, Y. Bian, L. Kang, P. L. Werner, and D. H. Werner, “Leap-frog continuous–discontinuous Galerkin time domain method for nanoarchitectures with the Drude model,” J. Lit. Technol. 35(22), 4888–4896 (2017). [CrossRef]

	Debye	Drude	Lorentz	Sellmeier	Critical Point
$a_{0}$	0	$ω_{D}^{2}$	$f_{L} ω_{L}^{2}$	$f_{S} ω_{S}^{2}$	$2 f_{c} ω_{c} (ω_{c} \cos φ_{c} - τ_{c} \sin φ_{c})$
$a_{1}$	$Δ ε τ^{- 1}$	0	0	0	$- 2 f_{c} ω_{c} \sin φ_{c}$
$b_{0}$	0	0	$ω_{L}^{2}$	$ω_{S}^{2}$	$ω_{c}^{2} + τ_{c}^{2}$
$b_{1}$	$τ^{- 1}$	$Γ_{D}$	$Γ_{L}$	0	$2 Γ_{c}$

	Number of Mesh Cells	Degrees of Freedom	CPU Time (s)	Memory (MB)
DGTD-GDM-PBCs	996	100,008	85	2,394
CST-FEM	34,234	226,454	230	1,320

	Debye	Drude	Lorentz	Sellmeier	Critical Point
$a_{0}$	0	$ω_{D}^{2}$	$f_{L} ω_{L}^{2}$	$f_{S} ω_{S}^{2}$	$2 f_{c} ω_{c} (ω_{c} \cos φ_{c} - τ_{c} \sin φ_{c})$
$a_{1}$	$Δ ε τ^{- 1}$	0	0	0	$- 2 f_{c} ω_{c} \sin φ_{c}$
$b_{0}$	0	0	$ω_{L}^{2}$	$ω_{S}^{2}$	$ω_{c}^{2} + τ_{c}^{2}$
$b_{1}$	$τ^{- 1}$	$Γ_{D}$	$Γ_{L}$	0	$2 Γ_{c}$

	Number of Mesh Cells	Degrees of Freedom	CPU Time (s)	Memory (MB)
DGTD-GDM-PBCs	996	100,008	85	2,394
CST-FEM	34,234	226,454	230	1,320

Discontinuous Galerkin time domain analysis of electromagnetic scattering from dispersive periodic nanostructures at oblique incidence

Abstract

1. Introduction

2. Methods and formulations

2.1. Transformed Maxwell’s equations and PBCs

2.2. Generalized dispersive material model

2.3. Discontinuous Galerkin spatial discretization

2.4. Runge-Kutta time stepping scheme

3. Numerical results

3.1. Reflection and transmission of a thin film

3.2. Surface plasmon excitation at the interface of a metal hole array

3.3. Absorption of a dual-band infrared absorber

4. Conclusions

Funding

References

Cited By

Figures (4)

Tables (2)

Equations (33)

Optics Express