Continuous-discontinuous Galerkin time domain (CDGTD) method with generalized dispersive material (GDM) model for computational photonics

Qiang Ren; Huaguang Bao; Sawyer D. Campbell; Ludmila J. Prokopeva; Alexander V. Kildishev; Douglas H. Werner

doi:10.1364/OE.26.029005

1. Introduction

The predictive power and computational efficiency of the numerical design tools for photonic applications heavily depend on the experimentally fitted dispersion models of the dielectric functions expressed in terms of a few physically meaningful, wavelength-independent parameters. While the Drude and Drude-Lorentz models remain the most popular physics-based dispersion approximations, a more efficient Drude-critical point model has also been introduced to fit the experimental data obtained from ellipsometric characterization of noble metals and alternative plasmonic materials [1–6]. Furthermore, various classical models have been proposed for other dispersive materials (e.g., dielectrics, semiconductor alloys), such as the Debye, Tauc-Lorentz, generalized oscillator, and generalized critical points [7] models along with some classical empirical models (e.g., the Cauchy and Cole-Cole models) [8]. In particular, an important class of dispersive materials in optics are plasmonic materials, in which an electric field excites the collective oscillations of free electrons (plasmons). Plasmonics has enabled new optical devices and architectures due to the ability of noble metals and alternative plasmonic materials to confine light to highly subwavelength scales [9,10]. Thus, due to the vast diversity of highly dispersive dielectric functions in nanophotonic applications, there is a need to have a robust universal approach that can model an arbitrary material dispersion ε(ω) in the time domain.

In contrast to frequency-domain numerical methods, where ε(ω) could be defined by a lookup table or a more general computable function, for which the Kramers-Krönig (KK) consistency is not necessarily enforced, time-domain methods require analytical causal representation in the time domain, e.g. in terms of a convolution integral or differential equations. This analytical description should allow effective discretization to calculate a local response at each point with an explicit recurrence. Such recurrences can be derived for the finite difference time domain (FDTD) and finite element time domain (FETD) methods using either a recursive convolution (RC) or an auxiliary differential equation (ADE) technique. Moreover, every dispersive model in the literature exploits its own type of RC or ADE approach, precluding formulation of a universal time integration scheme. This not only increases the complexity of the time-domain algorithm design, but also hinders the efficient maintenance of the code. Some representative work on time domain simulations of dispersive media can be found in [1–4,11–16]. To achieve a universal dispersion model for a wide range of dispersive media, a generalized dispersive material (GDM) model has recently been created based on [1/2] Padé approximants [14,15]. This model covers the Debye, Lorentz, Drude, Drude-Lorentz, and Drude-critical points (DCP) models with an arbitrary number of poles. The FDTD-GDM model has been applied to simulate light-matter interaction with gold nanostrips, gold nano-slits and multi-layer composite Ag/SiO2 films, semi-continuous metal films, graphene nanostructures, lasing in silver nanohole arrays, and its effectiveness has been validated against both analytical models and experimental data [14–19]. Also, a 2D version of the GDM model [19] was successfully implemented in DGTD simulations of graphene [20], where in contrast to this work, the dispersive surface conductivity of graphene is incorporated via numerical flux, while our method enables modeling of volumetric dispersive media without needing to modify the numerical flux.

In this paper, the volumetric GDM model has been incorporated, for the first time, into the CDGTD method [21–32] to enable numerical modeling of large, highly dispersive optical elements, including plasmonic materials. It inherits the advantages of the conventional DGTD method, such as an unstructured mesh, high-order elements, non-conformal meshing, and a nonspurious scheme based on the field variables E and B (EB scheme). Auxiliary differential equations are employed to couple the electric fields and polarization currents, while a leap-frog time integration scheme is proposed to solve Maxwell’s equations and ADEs with second-order accuracy. Finally, the integration of the GDM into the CDGTD method is validated against several analytical canonical problems. This work builds the basis for future research on applying DGTD/CDGTD-based methods to the efficient multiscale and multiphysics design of optical nanostructured materials.

2. Methods and formulations

2.1 Generalized dispersive material (GDM) model

The GDM model was first proposed in [14,16] using the [1/2] Padé approximants. The GDM model employed in this paper is a slight variation of the original formulation, in which the relationship between the polarization currents and the electric field can be cast into the following second-order ADE

\frac{\partial^{2} J_{m}}{\partial t^{2}} + b_{1, m} \frac{\partial J_{m}}{\partial t} + b_{0, m} J_{m} = a_{0, m} \frac{\partial E}{\partial t} + a_{1, m} \frac{\partial^{2} E}{\partial t^{2}}, m \in \bar{1, M},

where E is the electric field and

J_{m} = ε_{0}^{- 1} \partial P_{m} / \partial t

is the normalized polarization current from the m^th pole (or the m^th pair of poles). In the frequency domain the relationship between J and E (1) can be written as

J_{m} = \frac{- i ω}{ε_{0}} P_{m} = - i ω χ_{m} E,

where

χ_{m}

is the susceptibility in the frequency domain

χ_{m} (ω) = \frac{a_{0, m} - i ω a_{1, m}}{- ω^{2} - i ω b_{1, m} + b_{0, m}}

The full expression for the GDM relative permittivity

ε (ω)

consists of possibly multiple

χ_{m}

terms plus a finite constant permittivity

ε_{1}

and conductivity

σ

ε (ω) = ε_{1} - \frac{σ}{i ω ε_{0}} + \sum_{m} \frac{a_{0, m} - i ω a_{1, m}}{- ω^{2} - i ω b_{1, m} + b_{0, m}}

For example, the susceptibilities of the Debye, Drude, and Lorentz terms as well as the critical point pairs are

\begin{matrix} χ_{D e} (ω) = \frac{Δ ε / τ}{- i ω + 1 / τ} \\ χ_{D} (ω) = \frac{ω_{D}^{2}}{i ω (i ω - Γ_{D})}, \\ χ_{L} (ω) = \frac{f_{L} ω_{L}^{2}}{- ω^{2} - i ω Γ_{L} + ω_{L}^{2}} \\ χ_{c} (ω) = f_{c} ω_{c} (\frac{e^{i ϕ_{c}}}{ω_{c} - ω - i Γ_{c}} + \frac{e^{- i ϕ_{c}}}{ω_{c} + ω + i Γ_{c}}) = 2 f_{c} ω_{c} \frac{ω_{c} \cos ϕ_{c} + (i ω - Γ_{c}) \sin ϕ_{c}}{- ω^{2} - 2 i ω Γ_{c} + (ω_{c}^{2} + Γ_{c}^{2})} \end{matrix}

where Δε, τ, ω_D, Γ_D, f_L, ω_L, Γ_L, f_c, ω_c, Γ_c, and ϕ_c are constants whose physical meanings are described in [1,2,14]. For example, a large set of optical materials, specifically in plasmonics, can be characterized with a sum of the Drude term and the Lorentz (or critical point) terms. In the numerical implementation, all of the terms are expressed with universal parameters a_0,1, b_0,1 (3) that can be obtained from Table 1.

Table 1. Parameters of the ADE for Different Dispersion Terms

View Table | View all tables in this article

2.2 CDGTD method incorporating the GDM model

The first-order Maxwell’s equations with the polarization currents can be expressed as follows

ε_{1} \frac{\partial E}{\partial t} = c^{2} \nabla \times B - \frac{σ}{ε_{0}} E - J_{i} - \sum_{m = 1}^{M} J_{m},

\frac{\partial B}{\partial t} = - \nabla \times E - M_{i},

where E and B are the electric field intensity and the magnetic flux density, respectively; J_i is the imposed normalized electric current density, while M_i is the imposed magnetic current density; ε₀, c, ε₁, and σ denote the respective free-space permittivity, speed of light, constant relative permittivity, and finite electric conductivity (see (4)). Curl-conforming basis functions and divergence-conforming basis functions are employed to discretize the unknown field vectors E and B, respectively (EB scheme). For details on the construction process and the spatial distributions of these vector basis functions along with the advantages of the EB scheme, readers should refer to [21–27]. The polarization currents, J_m use the same curl-conforming basis functions and the same time grids as those for E. This configuration, which is referred to as EJ collocation, is widely used in both the DGTD and FDTD communities [2,27,33].

Next, a Galerkin scheme is used to test Maxwell’s Eqs. (6) and (7) with the GDM model (1), where curl-conforming basis functions Φ and divergence-conforming basis functions Ψ, are used, as in [21–23]. Additionally, an upwind flux scheme (Riemann solver) [22,23] is employed to decouple the adjacent subdomains, enabling the semi-discretized system of equations for the k^th subdomain to be written as

\begin{matrix} M_{e e}^{(k)} \frac{d e^{(k)}}{d t} = C_{e e}^{(k)} e^{(k)} + K_{e b}^{(k)} b^{(k)} + L_{e e}^{(k)} e^{(k)} + L_{e b}^{(k)} b^{(k)} + \sum_{l}^{} L_{e e}^{(k l)} e^{(l)}, \\ + \sum_{l}^{} L_{e b}^{(k l)} b^{(l)} + \sum_{m}^{} D_{m}^{(k)} j_{m}^{(k)} - j_{i}^{(k)} \end{matrix}

M_{b b}^{(k)} \frac{d b^{(k)}}{d t} = C_{b b}^{(k)} b^{(k)} + K_{b e}^{(k)} e^{(k)} + L_{b e}^{(k)} e^{(k)} + \sum_{l} L_{b b}^{(k)} b^{(l)} + \sum_{l} L_{b e}^{(k l)} e^{(l)} - m_{l}^{(k)},

\frac{\partial^{2} j_{m}^{(k)}}{\partial t^{2}} + b_{1, m} \frac{\partial j_{m}^{(k)}}{\partial t} + b_{0, m} j_{m}^{(k)} = a_{0, m} \frac{\partial e^{(k)}}{\partial t} + a_{1, m} \frac{\partial^{2} e^{(k)}}{\partial t^{2}},

where

D_{m, p, q}^{(k)} = ∭ ​ Φ_{​ p} \cdot Φ_{q}^{T} d V

, l is the index of the adjacent subdomain, V is the volume of the element, and T indicates a vector transpose. Definitions for the other system matrices can be found in [22,23].

Next, by applying central differences to (6) and (8) at $t_{n} = n Δ t$ and to (7) at $t_{n + 1 / 2} = (n + 1 / 2) Δ t$ we can obtain the system of linear equations as follows,

\begin{matrix} M_{b b}^{(k)} \frac{b_{n + 1 / 2}^{(k)} - b_{n - 1 / 2}^{(k)}}{Δ t} = C_{b b}^{(k)} \frac{b_{n + 1 / 2}^{(k)} + b_{n - 1 / 2}^{(k)}}{2} + K_{b e}^{(k)} e_{n}^{(k)} + L_{b b}^{(k)} b_{n - 1 / 2}^{(k)}, \\ + L_{b e}^{(k)} e_{n}^{(k)} + \sum_{l}^{} L_{b b}^{(k l)} b_{n - 1 / 2}^{(l)} + \sum_{l}^{} L_{b e}^{(k l)} e_{n}^{(l)} - m_{i, m}^{(k)} \end{matrix}

\begin{matrix} M_{e e}^{(k)} \frac{e_{n + 1}^{(k)} - e_{n}^{(k)}}{Δ t} = C_{e e}^{(k)} \frac{e_{n + 1}^{(k)} + e_{n}^{(k)}}{2} + K_{e b}^{(k)} b_{n - 1 / 2}^{(k)} + L_{e e}^{(k)} e_{n}^{(k)} + L_{e b}^{(k)} b_{n + 1 / 2}^{(k)} \\ + \sum_{l}^{} L_{e e}^{(k l)} e_{n}^{(l)} + \sum_{l}^{} L_{e b}^{(k l)} b_{n + 1 / 2}^{(l)} + \sum_{m}^{} D_{m}^{(k)} \frac{j_{m, n + 1}^{(k)} + j_{m, n}^{(k)}}{2} - j_{i, n + 1 / 2}^{(k)}, \end{matrix}

\begin{matrix} \frac{j_{m, n + 1}^{(k)} - 2 j_{m, n}^{(k)} + j_{m, n - 1}^{(k)}}{{(Δ t)}^{2}} + b_{1} \frac{j_{m, n + 1}^{(k)} - j_{m, n - 1}^{(k)}}{2 Δ t} + b_{0} j_{m, n}^{(k)} = \\ a_{0} \frac{e_{n + 1}^{(k)} - e_{n - 1}^{(k)}}{2 Δ t} + a_{1} \frac{e_{n + 1}^{(k)} - 2 e_{n}^{(k)} + e_{n - 1}^{(k)}}{{(Δ t)}^{2}} \end{matrix}

After merging similar terms in (13),

j_{m, n + 1}^{(k)}

can be expressed as

j_{m, n + 1}^{(k)} = α_{m} j_{m, n - 1}^{(k)} + β_{m} j_{m, n - 1}^{(k)} + γ_{m} e_{n + 1}^{(k)} + ξ_{m} e_{n}^{(k)} + ζ_{m} e_{n - 1}^{(k)},

where the coefficients α_m, β_m, γ_m, ξ_m, ζ_m are given by

\begin{array}{l} α_{m} = {[1 + \frac{b_{1} Δ t}{2}]}^{- 1} [2 - b_{0} {(Δ t)}^{2}], β_{m} = {[1 + \frac{b_{1} Δ t}{2}]}^{- 1} [\frac{b_{1} Δ t}{2} - 1], \\ γ_{m} = {[1 + \frac{b_{1} Δ t}{2}]}^{- 1} [a_{1} + \frac{a_{0} Δ t}{2}], ξ_{m} = {[1 + \frac{b_{1} Δ t}{2}]}^{- 1} (- 2 a_{1}), \\ ζ_{m} = {[1 + \frac{b_{1} Δ t}{2}]}^{- 1} [a_{1} - \frac{a_{0} Δ t}{2}] . \end{array}

The term

b_{n + 1 / 2}^{(k)}

in (11) can be updated as

\begin{matrix} b_{n + 1 / 2}^{(k)} = {[M_{b b}^{(k)} - \frac{Δ t}{2} C_{b b}^{(k)}]}^{- 1} [(M_{b b}^{(k)} + \frac{Δ t}{2} C_{b b}^{(k)}) b_{n - 1 / 2}^{(k)} \\ + Δ t (K_{b e}^{(k)} e_{n}^{(k)} + L_{b b}^{(k)} b_{n - 1 / 2}^{(k)} + L_{b e}^{(k)} e_{n}^{(k)} + \sum_{l}^{} L_{b b}^{(k l)} b_{n - 1 / 2}^{(l)} + \sum_{l}^{} L_{b e}^{(k l)} e_{n}^{(l)} - m_{i, n}^{(k)})] . \end{matrix}

Substituting (14) into (12), it follows that

e_{n + 1}^{(k)}

can be obtained via

\begin{matrix} e_{n + 1}^{(k)} = {[M_{e e}^{(k)} \frac{Δ t C_{e e}^{(k)}}{2} - \frac{Δ t \sum_{m = 1}^{M} γ_{m} D_{m}^{(k)} Z_{s l} Z_{l s}}{2}]}^{- 1} \times \\ {[M_{e e}^{(k)} \frac{Δ t C_{e e}^{(k)}}{2} - \frac{Δ t \sum_{m = 1}^{M} ξ_{m} D_{m}^{(k)} Z_{s l} Z_{l s}}{2}]} e_{n}^{(k)} \\ + Δ t (K_{e b}^{(k)} b_{n - 1 / 2}^{(k)} + L_{e e}^{(k)} e_{n}^{(k)} + L_{e b}^{(k)} b_{n + 1 / 2}^{(k)} + \sum_{l} L_{b b}^{(k l)} b_{n - 1 / 2}^{(k)} + \sum_{l} L_{b e}^{(k l)} e_{n}^{(l)} - m_{i, n}^{(k)}) \\ - Δ t \sum_{m = 1}^{M} D_{m}^{(k)} (\frac{ξ_{m}}{2} Z_{s l} Z_{l s} e_{n - 1}^{(k)} + \frac{1 + α_{m}}{2} Z_{s l} J_{m, n}^{(k)} + \frac{β_{m}}{2} Z_{s l} J_{m, n - 1}^{(k)})} . \end{matrix}

Then

j_{m, n + 1}^{(k)}

can be updated according to

j_{m, n + 1}^{(k)} = α_{m} j_{m, n}^{(k)} + β_{m} j_{m, n - 1}^{(k)} + Z_{l s} (γ_{m} e_{m, n + 1}^{(k)} + ξ_{m} e_{m, n}^{(k)} + ζ_{m} e_{m, n - 1}^{(k)}) .

Since the dispersive media may not occupy the entire volume of the k^th subdomain, the unknown vectors

e^{(k)}

and

j_{m}^{(k)}

may have different lengths. Thus, the matrix

Z_{s l}

,

(Z_{l s} = Z_{s l}^{⊺})

is required to map the indices of unknowns for the polarization currents to the DoF indices of the k^th subdomain.

3. Numerical results

Conventionally, the dispersion of metals is described with Drude-Lorentz models, where models employing up to ten Lorentz terms have been presented in the literature. However, recent developments have yielded more accurate fits with just two oscillators by adding a non-zero phase shift to the damped Lorentz oscillator. In the numerical tests, we use gold, as an example of a highly dispersive plasmonic material. Figure 1(b) shows a depiction of gold’s complex permittivity along with experimentally measured data [34]. The solid lines are the results calculated with the Drude term plus a two critical points (D2CP) dispersive model and the dashed lines indicate the experimental data. Here, for demonstration of the proposed numerical approach we show broadband simulation results of a gold film and scattering from a gold sphere. The dispersion parameters of gold used in the simulations are summarized in Table 2, where the oscillation frequencies and damping parameters are shown in nanometers with $2 π c / λ_{m}$ and $2 π c / Γ_{m}$ being the conversion formulas to their original values in rad/s, respectively.

Fig. 1 Simulation setup of a D2CP film under plane wave illumination. (a) The geometry, source, and probes for the computational domain configuration. (b) Permittivity of gold with D2CP model compared to experimental data [34]. (c) The transient incident wave. (d,e) The transient reflected (d) and transmitted (e) E field at probes 1 and 2, respectively.

Download Full Size | PDF

Table 2. The Parameters of Drude - 2 Critical Points Model

View Table | View all tables in this article

3.1 Reflection and transmission of a thin film

The configuration of the computational domain for this test case is illustrated in Fig. 1(a). The 50-nanometer-thick gold film is described with a D2CP dispersive model [14] and listed in Table 2. The incident field is generated by a modulated Gaussian pulse

E_{pulse} (t) = e^{- {(\frac{t - t_{0}}{t_{w}})}^{2}} \sin [ω (t - t_{0})] .

where the carrier wavelength,

λ = 2 π c / ω

, is at 366 nm, the pulse width t_w, used to probe the response numerically in the optical range, from 200 nm to 2000 nm, is 0.77 fs, and the pulse center is

t_{0} = 3 t_{w}

. This source is located at x_s = −2.4 μm and thus generates the incident plane wave

E^{inc} (x, t) = E_{pulse} (t - (x - x_{s}) / c)

propagating along the x direction, and the electric field is polarized along the y direction.

For the detailed implementation of the plane wave excitation for the vector basis functions, please refer to [26]. In the x direction, two perfectly matched layers (PMLs) are used to truncate the computational domain. In the lateral direction, the boundaries perpendicular to the y and z axes are assigned as PEC and PMC to emulate periodic boundary conditions. The entire region is meshed with hexahedral elements and second order spatial basis functions were employed. Here, the term points per film (ppf) is introduced as a measure to quantify the mesh size, which is equal to 50/ppf nanometers. The time step used for the simulations is determined by $Δ t = 0.16 Δ x / c$ , where a small Courant number is taken because of stability limitations of the discontinuous Galerkin scheme. It should be noted that the incorporation of the GDM model has no impact on the stability of the conventional CDGTD method. Finally, two probes are placed at the indicated locations to record the transient reflected and transmitted electric field.

The relative error of the reflection and transmission coefficients is used to test the accuracy of the GDM-CDGTD implementation, which is defined as

Error of X = \frac{| X_{n} (λ) - X_{a} (λ) |}{| X_{a} (λ) |} \times 100 %,

where X is the reflection or transmission coefficient, the subscript n represents the numerical result, and a denotes the analytical solution. First, the electric fields at the probes plotted in Fig. 1(d,e) are recorded for a ppf value of 2. The reflection and transmission coefficients are calculated using FFT over the wavelength range of [200 nm, 2000 nm] and compared to the analytical solution, Fig. 1(a,b). Second, the convergence of the proposed method is studied. We perform simulations for different mesh sizes, in which the element sizes are 50 nm (ppf = 1), 25 nm (ppf = 2), 12.5 nm (ppf = 4), and 6.25 nm (ppf = 8), respectively.

Figure 2(c), 2(d) shows the convergence of the relative errors of the reflection and transmission coefficients with respect to mesh size. As a reference, we also plot numerical errors using a similar setup for the standard Yee’s FDTD and a commercial FETD solver with linear finite elements (both with ppf = 8) employing the same GDM model as in the CDGTD method. To check the order of convergence, in Fig. 2(e), 2(f) we then plot the ratio of the relative error for ppf/2 to that for a given ppf value. Again, as a reference we plot the convergence ratio for the standard second order FDTD and FETD methods, which just slightly deviate numerically from a theoretical value of 4 for all ppf values (the plot shows a case of ppf = 8).

Fig. 2 The reflection and transmission coefficients and their relative errors for different ppf values. (a), (b) Comparison of the reflection and transmission coefficients between the analytical solution and the CDGTD method with a GDM model for gold and ppf = 2. (c), (d) Relative errors of the reflection and transmission coefficients. (e), (f) Ratio of the relative error for ppf/2 to that for ppf.

Download Full Size | PDF

The analysis indicated the following advantages of the proposed method, which are essential for efficient design and optimization with direct solvers. First, the results in Fig. 2(a), 2(b) show that just 2 cells per film are already sufficient for a 1% accurate simulation over the entire range of wavelengths, while reference methods conventionally used in computational nanophotonics require at least a 4 times finer mesh, Fig. 2(c), 2(d). For meshes with ppf = 2 and 4, the proposed method exhibits convergence behavior beyond the second-order methods (FDTD and FETD), and then it saturates for extremely refined meshes (ppf = 8). This saturated behavior is a consequence of the spatial basis functions employed in the simulation and is consistent with previous studies reported in the literature [35]. While these basis functions may lead to a potential increase in error when the mesh density becomes too fine, they also enable elements with lengths of up to two wavelengths in contrast to traditional low-order basis functions, which typically limit element lengths to 0.1 wavelengths.

Therefore, the proposed method is advantageous versus the conventional techniques for topological optimization problems, where massive and diverse shape changes are spanned. In this case, its ability to perform accurate (less than 1% error) simulations on a coarse mesh can result in extreme improvement in computational efficiency.

3.2 Near-field scattering from a nanosphere

The near-field response of an incident plane wave scattering by a gold nanosphere, whose material parameters are also listed in Table 2, is evaluated with the proposed method to provide a further demonstration of its validity. The radius of the nanosphere is 150 nm. The plane wave propagation direction vector is [0, 0, 1], and the electric field polarization vector is [1, 0, 0]. The source type of the wideband excitation plane wave is a 1st order Blackman-Harris window (BHW) pulse [36] with a characteristic frequency of 400 THz (wavelength is 750 nm). The total fields are recorded at three different observation points located at [-250, −250, −250] nm, [-250, −250, 0] nm and [-250, −250, 250] nm, respectively, as illustrated in Fig. 3(a). The physical region is divided into 20 subdomains in total. Here, the proposed CDGTD method is based on $E$ and $B$ fields, instead of $E$ and $H$ and, consequently, the nanosphere and its surrounding space are meshed with tetrahedrons employing linear-tangential/quadratic-normal (LT/QN) basis functions for the electric field intensity and linear-normal/quadratic-tangential (LN/QT) basis functions for the magnetic flux density [23] as shown in Fig. 3(b). The outer region and the perfectly matched layer (PML) are meshed with high-order brick elements (not shown).

Fig. 3 Geometry and mesh of the nanosphere. (a) Geometry of the nanosphere with source and observation points. (b) Tetrahedral mesh of the nanosphere and its surrounding space.

Download Full Size | PDF

The total electric fields recorded at the observation points are transformed to the frequency domain and normalized to the incident fields. Figure 4 demonstrates good agreement between the calculated results from the CDGTD method employing GDM and the analytical solutions obtained from the Mie theory [37,38] for the frequency range of [300 THz, 900 THz] (wavelength range of [333.3 nm, 1000 nm]), which also confirms the effectiveness of the proposed method.

Fig. 4 Near-field response of the gold nanosphere employing a DCP dispersive model. The numerical results from the CDGTD method with the GDM model are compared to the analytical solutions from Mie theory. (a) The total fields at [-250, −250, −250] nm. (b) The total fields at [-250, −250, 0] nm. (c) The total fields at [-250, −250, 250] nm.

Download Full Size | PDF

4. Conclusions

In this paper, a generalized dispersive material model has been successfully incorporated into the CDGTD method using an ADE approach. This multiphysics coupling enables the simulation of any linear optically dispersive media for nanophotonic applications, including, for example, plasmonic and polaritonic materials, where accurate broadband numerical modeling of the material dispersion in the time domain is essential. The advantages of the proposed method versus conventional approaches are shown with reference tests by demonstrating that CDGTD with an embedded GDM model provides accurate simulations on 4 times coarser meshes. The technique is implemented for a general 3D formulation and is also tested for a classic case of scattering from a dispersive metal sphere. Coupling the GDM model with an accurate, easily scalable CDGTD computational electromagnetics engine makes it applicable to a vast range of design and optimization problems in nanophotonics. Our study opens up a pathway for future applications of computationally efficient DGTD/CDGTD-based techniques to the multiscale and multiphysics design and optimization of novel nanostructured photonic elements and optoelectronic systems.

Funding

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

References

1. 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 Propagat. 60(4), 1969–1977 (2012). [CrossRef]

2. A. Taflove and S. C. Hagness, Computational Electromagnetics: The Finite-Difference Time-Domain Method (Artech House, 2000).

3. F. Hao and P. Nordlander, “Efficient dielectric function for FDTD simulation of the optical properties of silver and gold nanoparticles,” Chem. Phys. Lett. 446(1-3), 115–118 (2007). [CrossRef]

4. N. Theethayi, Y. Baba, F. Rachidi, and R. Thottappillil, “On the choice between transmission line equations and full-wave Maxwell’s equations for transient analysis of buried wires,” IEEE Trans. Electromagn. Compat. 50(2), 347–357 (2008). [CrossRef]

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

6. P. G. Etchegoin, E. C. L. Ru, and M. Meyer, “Erratum: An analytic model for the optical properties of gold,” J. Chem. Phys. 127(18), 189901 (2007). [CrossRef]

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

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

9. M. I. Stockman, K. Kneipp, S. I. Bozhevolnyi, S. Saha, A. Dutta, J. Ndukaife, N. Kinsey, H. Reddy, U. Guler, V. M. Shalaev, A. Boltasseva, B. Gholipour, H. N. S. Krishnamoorthy, K. F. MacDonald, C. Soci, N. I. Zheludev, V. Savinov, R. Singh, P. Groß, C. Lienau, M. Vadai, M. L. Solomon, D. R. Barton III, M. Lawrence, J. A. Dionne, S. V. Boriskina, R. Esteban, J. Aizpurua, X. Zhang, S. Yang, D. Wang, W. Wang, T. W. Odom, N. Accanto, P. M. de Roque, I. M. Hancu, L. Piatkowski, N. F. van Hulst, and M. F. Kling, “Roadmap on plasmonics,” J. Opt. 20(4), 043001 (2018). [CrossRef]

10. S. A. Maier, Plasmonics: Fundamentals and Applications (Springer Science & Business Media, 2007).

11. J. Y. Lu and Y. H. Chang, “Optical singularities associated with the energy flow of two closely spaced core-shell nanocylinders,” Opt. Express 17(22), 19451–19458 (2009). [CrossRef] [PubMed]

12. A. Vial and T. Laroche, “Description of dispersion properties of metals by means of the critical points model and application to the study of resonant structures using the FDTD method,” J. Phys. D Appl. Phys. 40(22), 7152–7158 (2007). [CrossRef]

13. K. P. Prokopidis and D. C. Zografopoulos, “A unified FDTD/PML scheme based on critical points for accurate studies of plasmonic structures,” J. Lightwave Technol. 31(15), 2467–2476 (2013). [CrossRef]

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

15. L. J. Prokopeva, A. V. Kildishev, J. Fang, J. Borneman, M. D. Thoreson, V. M. Shalaev, and V. P. Drachev, “Nanoplasmonics FDTD simulations using a generalized dispersive material model,” Proceedings of 27th International Review of Progress in Applied Computational Electromagnetics Symposium (ACES 2011), pp. 1–6.

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

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

18. S. I. Azzam, J. Fang, J. Liu, Z. Wang, N. Arnold, T. Klar, L. J. Prokopeva, X. Meng, V. M. Shalaev, and A. V. Kildishev, “Time-resolved dynamics of plasmonic systems with gain,” (unpublished).

19. L. J. Prokopeva and A. V. Kildishev, “Efficient time-domain model of the graphene dielectric function,” in [SPIE NanoScience+ Engineering], 88060H, International Society for Optics and Photonics (2013).

20. J. F. M. Werra, C. Wolff, C. Matyssek, and K. Busch, “Current sheets in the discontinuous Galerkin time-domain method: an application to graphene,” Proc. SPIE 9502, 95020E (2015). [CrossRef]

21. Q. Ren, Q. Sun, L. Tobón, Q. Zhan, and Q. H. Liu, “EB scheme based hybrid SE-FE DGTD method for multiscale EM simulations,” IEEE. Trans. Antennas Propagat. 64(9), 4088–4091 (2016). [CrossRef]

22. Q. Ren, L. E. Tobón, Q. Sun, and Q. H. Liu, “A new 3D non-spurious discontinuous Galerkin spectral element time domain (DG-SETD) method for Maxwell’s equations,” IEEE. Trans. Antennas Propagat. 63(6), 2585–2594 (2015). [CrossRef]

23. L. E. Tobón, Q. Ren, and Q. H. Liu, “A new efficient 3D discontinuous Galerkin time domain (DGTD) method for large and multiscale electromagnetic simulations,” J. Comput. Phys. 283, 374–387 (2015). [CrossRef]

24. Q. Ren, L. E. Tobón, and Q. H. Liu, “A new 2D non-spurious discontinuous Galerkin finite element time domain (DG-FETD) method for Maxwell’s equations,” Prog. Electromagnetics Res. 143, 385–404 (2013). [CrossRef]

25. Q. Sun, Q. Zhan, Q. Ren, and Q. H. Liu, “Wave equation-based implicit subdomain DGTD method for modeling of electrically small problems,” IEEE Trans. Microw. Theory Tech. 65(4), 1111–1119 (2017). [CrossRef]

26. Q. Ren, Q. Zhan, and Q. H. Liu, “An improved subdomain level non-conformal discontinuous Galerkin time domain (DGTD) method for materials with full-tensor constitutive parameters,” IEEE Photonics J. 9(2), 2600113 (2017). [CrossRef]

27. 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. Lightwave Technol. 35(22), 4888–4896 (2017). [CrossRef]

28. J. Chen, L. E. Tobón, M. Chai, J. A. Mix, and Q. H. Liu, “Efficient implicit-explicit time stepping scheme with domain decomposition for multiscale modeling of layered structures,” IEEE Trans. Compon. Packaging Manuf. Technol. 1(9), 1438–1446 (2011). [CrossRef]

29. L. D. Angulo, J. Alvarez, F. L. Teixeira, M. F. Pantoja, and S. G. Garcia, “A nodal continuous-discontinuous Galerkin time-domain method for Maxwell’s equations,” IEEE Trans. Microw. Theory Tech. 63(10), 3081–3093 (2015). [CrossRef]

30. J. Chen, Q. H. Liu, M. Chai, and J. A. Mix, “A non-spurious 3-D vector discontinuous Galerkin finite-element time-domain method,” IEEE Microw. Wirel. Compon. Lett. 20(1), 1–3 (2010). [CrossRef]

31. S. Gedney, C. Luo, J. Roden, R. Crawford, B. Guernsey, J. Miller, and E. Lucas, “A discontinuous Galerkin finite element time domain method with PML,” Proceedings of IEEE Antennas and Propagation Society International Symposium (IEEE 2008), pp. 1–4.

32. F. G. Hu and C. F. Wang, “Modeling of waveguide structures using DG-FETD method with higher-order tetrahedral elements,” IEEE Trans. Microw. Theory Tech. 60(7), 2046–2054 (2012). [CrossRef]

33. Y. Yu and J. J. Simpson, “An E-J collocated 3-D FDTD model of electromagnetic wave propagation in magnetized cold plasma,” IEEE Trans. Antenn. Propag. 58(2), 469–478 (2010). [CrossRef]

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

35. M. M. Ilic and B. N. Notaros, “Higher order hierarchical curved hexahedral vector finite elements for electromagnetic modeling,” IEEE Trans. Microw. Theory Tech. 51(3), 1026–1033 (2003). [CrossRef]

36. Q. H. Liu, “An FDTD algorithm with perfectly matched layers for conductive media,” Microw. Opt. Technol. Lett. 14(2), 134–137 (1997). [CrossRef]

37. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley and Sons, 2008).

38. J. Schafer, S. Lee, and A. Kienle, “Calculation of the near fields for the scattering of electromagnetic waves by multiple infinite cylinders at perpendicular incidence,” J. Quant. Spectrosc. Radiat. Transf. 113(16), 2113–2123 (2012). [CrossRef]

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

	ε₁	f_m	λ_m [nm]	Γ_m [nm]	φ_m [deg]
D	2.75	N/A	140	17000	N/A
C₁		0.75	445	2400	−22
C₂		1.55	140	950	−38.6

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

	ε₁	f_m	λ_m [nm]	Γ_m [nm]	φ_m [deg]
D	2.75	N/A	140	17000	N/A
C₁		0.75	445	2400	−22
C₂		1.55	140	950	−38.6

Continuous-discontinuous Galerkin time domain (CDGTD) method with generalized dispersive material (GDM) model for computational photonics

Abstract

1. Introduction

2. Methods and formulations

2.1 Generalized dispersive material (GDM) model

2.2 CDGTD method incorporating the GDM model

3. Numerical results

3.1 Reflection and transmission of a thin film

3.2 Near-field scattering from a nanosphere

4. Conclusions

Funding

References

Cited By

Figures (4)

Tables (2)

Equations (20)

Optics Express