Conductive mixed-order generalized dispersion model for noble metals in the optical regime

Wending Mai; Sawyer D. Campbell; Douglas H. Werner

doi:10.1364/OE.435297

1. Introduction

Metals are highly dispersive in the optical regime and numerous experiments have been performed to measure their material properties at different frequencies. In order to simulate dispersive materials in time-domain solvers, such frequency dependent properties have to be expressed in closed form models [1]. Some of these dispersive models have a direct physical meaning. For example, the Debye model captures relaxation [2–4], the Lorentz model is used for resonance [4], the Sellmeier model is used to represent narrow-band absorption [5], and the Drude model is characteristic of cold-plasma performance [1,4]. On the other hand, many dispersive models are proposed without clear physical meaning, such as the Critical Points model [6–8]. Although devoid of physical explanation, those mathematical expressions have their own merits with more freedom for accurate fitting. Moreover, modern developments of artificial metamaterials have engineered dispersive behaviors beyond conventional physical interpretations, and thus can only be represented by mathematical dispersion models. To accurately fit various materials in a unified approach, the Generalized Dispersion Model (GDM) was proposed as the sum of a series of Padé polynomials [9,10]. The GDM encompasses a large number of different dispersive models, both those that are fundamentally rooted in physics and others that are purely mathematical in nature. Therefore, it possess the ability to model almost arbitrary material dispersion behavior [9–13].

Regardless of their specific formulation, these closed-form dispersive models are implemented into time-domain solvers via convolution or through an auxiliary differential equation (ADE) method [1]. The latter is simpler to derive and can be easily generalized to nonlinear cases. To build up the ADE scheme for a time-domain solver, the time derivative terms are replaced by auxiliary vectors. Consequently, higher order derivative equations require additional vectors which increase CPU and memory consumption burdens [4].

Researchers have found the Drude-2-critical points (D2CP) model to be superior to the Drude and Drude-Lorentz models for approximating the dispersion of metals in the optical regime and therefore it is widely used in time domain solvers [9,10,14]. The D2CP model contains a Drude term and two critical point terms and has been incorporated into conventional GDM, which is composed of three second-order polynomials. Consequently, there are 8 unknown vectors ($\vec{E}$, $\vec{H}$, and 6 auxiliary vectors) associated with the conventional GDM D2CP model.

In this paper, we propose a conductive mixed-order general dispersive model (CMGDM). Similar to the conventional GDM model, in CMGDM different materials share the same model expression albeit with different parameter sets. We found that by adding an extra conductive term and mixing with a 1^st order Padé polynomial, the differential order of the Drude term can also be reduced. The proposed CMGDM is mathematically equivalent to the traditional GDM and therefore has identical accuracy. However, for materials containing a Drude term, CMGDM has fewer unknowns when implemented with an ADE method. In practice, most metals in the optical regime can be accurately represented by a D2CP model that contains a 1^st order Drude term. Therefore, efficiency improvements can be expected for most metallo-dielectric device simulations in the optical regime. Specifically, compared to the conventional GDM with only a Debye or Drude term, the total unknown vectors of CMGDM are reduced from 4 to 3, representing a reduction of about 25%. Even for a more complex D2CP model, the total number of unknown vectors in a CMGDM model is reduced from 8 to 7. Therefore, the computational time and memory consumption can be reduced by about 12.5% compared to a conventional GDM implementation.

In the following sections, we first revisit the conventional GDM D2CP method and derive the GDM parameters of several noble metals: Au, Ag, and Al. Then we propose the CMGDM model with a conductive term and mixed-order Padé polynomials. Finally, CMGDM models are integrated with a DGTD solver to validate their performance improvements over the conventional GDM model.

2. Conventional D2CP GDM model

2.1 Parameters of GDM model

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

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

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

where the component $j\omega $ in the frequency domain represents the engineering convention corresponding to time-varying fields as ${e^{j\omega }}$. The dispersive permittivity can be modeled by a conventional GDM model composed of the permittivity ${\varepsilon _r}$ at the infinite high frequency limit (${\varepsilon _\infty })$ and M partial terms [20]:

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

where ${\varepsilon _\infty }$ is a non-dispersive term while the second term represents the frequency-dependent behavior. The dispersive term here is the sum of M 2^nd order Padé polynomials. Its dispersion can be rearranged by manipulating the parameters ${a_0}$, ${a_1}$, $\textrm{\; }{b_0}$, and $\textrm{\; }{b_1}$ to represent multiple classical dispersion formulas such as the Debye (De), Drude (Dr), Lorentz (L), Sellmeier (S), and critical points (C) models whose dispersive terms are expressed as [9,15]:

(4)$${\varepsilon _{De}}(\omega )= \frac{{{C_{De}}}}{{1 + j\omega {\tau _{De}}}}$$

(5)$${\varepsilon _{Dr}}(\omega )={-} \frac{{\omega _{Dr}^2}}{{\omega ({\omega - j{\tau_{Dr}}} )}}$$

(6)$${\varepsilon _L}(\omega )= \frac{{{C_L}\omega _L^2}}{{\omega _L^2 - {\omega ^2} + j\omega {\tau _L}}}$$

(7)$${\varepsilon _S}(\omega )= \frac{{{C_S}\omega _S^2}}{{\omega _S^2 - {\omega ^2}}}$$

(8)$${\varepsilon _C}(\omega )= {C_C}\left( {\frac{{\cos (\varphi )+ jsin(\varphi )}}{{{\omega_C} - \omega - j{\tau_C}}} + \frac{{\cos (\varphi )- jsin(\varphi )}}{{{\omega_C} + \omega + j{\tau_C}}}} \right),$$

where ${\omega _{Dr}}$, ${\omega _L}$, ${\omega _S}$, $\; {\omega _c}$, ${\tau _{De}}$, ${\tau _{Dr}}$, ${\tau _L}$, ${\tau _C}$, ${C_{De}}$, ${C_L}$, ${C_S}$ and ${C_C}$ are constant parameters. By comparing Eqs. (4–8) with Eq. (3), a relationship between these parameters and the conventional GDM parameters ${a_0}$, ${a_1}$, ${b_0}$, and ${b_1}$ can be found (see Table 1):

The Drude and Critical Points terms in Eq. (5) and Eq. (8) are given in their frequency dependent form. Sometimes, it is convenient also to change to a wavelength dependent dielectric function which is more commonly utilized in optics and spectroscopy.

(9)$${\varepsilon _{Dr}}(\lambda )={-} \frac{1}{{\lambda _{Dr}^2\left( {\frac{1}{{{\lambda^2}}} - \frac{j}{{{\gamma_{Dr}}\lambda }}} \right)}}, $$

(10)$${\varepsilon _C}(\lambda )= \frac{A}{{{\lambda _c}}}\left( {\frac{{\cos (\varphi )+ jsin(\varphi )}}{{\frac{1}{{{\lambda_c}}} - \frac{1}{\lambda } - \frac{j}{{{\gamma_C}}}}} + \frac{{\cos (\varphi )- jsin(\varphi )}}{{\frac{1}{{{\lambda_c}}} + \frac{1}{\lambda } + \frac{j}{{{\gamma_C}}}}}} \right), $$

where $\lambda = 2\pi c/\omega $ (with c = speed of light), ${\lambda _{Dr}} = 2\pi c/{\omega _{Dr}}$ (plasma wavelength), ${\gamma _{Dr}} = 2\pi c/{\tau _{Dr}}$ (damping expressed in wavelengths), ${\lambda _C} = 2\pi c/{\omega _C}$ (interband transition wavelength), ${\gamma _C} = 2\pi c/{\tau _C}$ (transition broadening expressed in wavelengths), and $A = {C_c}/\omega $ (dimensionless critical points amplitude).

Table 1. Conventional GDM parameters for various dispersion terms^a

View Table | View all tables in this article

In order to model metals in the violet/near-UV region, interband transitions have to be taken into consideration. For example, Au has at least two interband transitions at $\lambda \sim 470$ and ${\sim} 330\; nm$. Their line shapes are not very well accounted for by a Lorentz oscillator. Therefore, adding a Lorentz term to a Drude term faces accuracy limitations when modeling interband transitions. On the contrary, the Lorentzian-Gaussian model has higher accuracy for interband transitions as it naturally satisfies the Kramers-Kronig relations [16]. As a compromise, the critical points term can achieve better accuracy than the Drude-Lorentz model, while on the other hand, it is convertible to a GDM model with better flexibility than the Lorentzian-Gaussian model. Therefore, an analytical model with a Drude term and two Critical Points terms is widely used to represent wideband dispersion of metals [6]. Researchers have also employed D2CP to model other dispersive metals, including Ag, Al, Cu, Cr, etc. [6,7,14,17,18]. The corresponding parameters for GDM can be deduced through the above derivation. Further investigation found that the D2CP model given for Cu and Cr works only over a narrow band [7]. So here we only list recommended GDM parameters for Au [6], Ag [7], and Al [7] for comparatively wide band accuracy (Table 2).

Table 2. Recommended conventional GDM parameters for noble metals

View Table | View all tables in this article

2.2 ADE implement of conventional GDM model

In order to implement the conventional D2CP GDM into a time domain solver, substituting Eq. (3) into Eq. (2) we have:

(11)$$\nabla \times \vec{H} = \sigma \vec{E} + {\varepsilon _0}{\varepsilon _\infty }j\omega \vec{E} + {\varepsilon _0}j\omega \mathop \sum \nolimits_m {\vec{P}_m}$$

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

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

We note here that Eq. (12) has an ${\omega ^2}$ term. Hence, in order to solve it by using the ADE method, an additional auxiliary term ${\vec{Q}_m}$ was introduced to lower the order [11–13], where

(13)$${\vec{Q}_m} = j\omega {\vec{P}_m}. $$

Equation (12) can therefore be rewritten as:

(14)$$j\omega {\vec{Q}_m} + {b_{1,m}}{\vec{Q}_m} + {b_{0,m}}{\vec{P}_m} = {a_{0,m}}\vec{E} + {a_{1,m}}j\omega \vec{E}. $$

As shown, to solve the Maxwell equations based on a conventional D2CP GDM model, 8 unknown vectors have to be solved ($\vec{E}$, $\vec{H}$, ${\vec{P}_1}$, ${\vec{P}_2}$, ${\vec{P}_3}$, ${\vec{Q}_1}$, ${\vec{Q}_2}$, ${\vec{Q}_3}$). Among them, ${\vec{P}_1}$ and ${\vec{Q}_1}$ correspond to the Drude term, while ${\vec{P}_2}$, ${\vec{Q}_2}$ and ${\vec{P}_3}$, ${\vec{Q}_3}$ are associated with the 1^st and 2^nd critical points terms.

3. Conductive mixed-order GDM model

3.1 Derivation

Instead of fitting the dispersion via a 2^nd order Padé polynomial, here we employ a model with mixed order. Suppose we have N 1^st order terms, and M 2^nd order terms, such that

(15)$${\varepsilon _r}(\omega )= {\varepsilon _\infty } + \mathop \sum \nolimits_{n = 1}^N \frac{{{a_{0,n}}}}{{{b_{0,n}} + j\omega }} + \mathop \sum \nolimits_{m = N + 1}^M \frac{{{a_{0,m}} + j\omega {a_{1,m}}}}{{{b_{0,m}} + j\omega {b_{1,m}} - {\omega ^2}}}, $$

where ${a_{0,n}}$ and ${b_{0,n}}$ are employed to represent 1^st order dispersion formulas.

Obviously, the 1^st order Padé polynomial will reduce to the Debye model in Eq. (4) by setting:

(16)$$\left\{ {\begin{array}{*{20}{c}} {{a_{0,n}} = \frac{{{C_{De}}}}{{{\tau _{De}}}}}\\ {{b_{0,n}} = \frac{1}{{{\tau _{De}}}}} \end{array}} \right.$$

On the other hand, lowering the order of the Drude model is comparatively more complicated. To this end, we introduced an extra term $\frac{C}{{j\omega }}$ where C is a constant. By summing up the 1^st order Padé polynomial and the extra term, we have:

(17)$$\frac{{{a_{0,n}}}}{{{b_{0,n}} + j\omega }} + \frac{{{C_n}}}{{j\omega }} = \frac{{{a_{0,n}}j\omega + C{b_{0,n}} + Cj\omega }}{{\omega ({\omega - {b_{0,n}}j} )}}. $$

Next, by comparing with Eq. (5), we arrive at:

(18)$${\varepsilon _{Dr}}(\omega )= \frac{{{a_{0,n}}}}{{{b_{0,n}} + j\omega }} - \frac{{{a_{0,n}}}}{{j\omega }} = \frac{{{a_{0,n}}{b_{0,n}}}}{{\omega ({\omega - {b_{0,n}}j} )}} ={-} \frac{{\omega _{Dr}^2}}{{\omega ({\omega - j{\tau_{Dr}}} )}}, $$

where

(19)$$\left\{ {\begin{array}{{c}} {{a_{0,n}} ={-} \frac{{\omega_{Dr}^2}}{{{\tau_{Dr}}}}}\\ {{b_{0,n}} = {\tau_{Dr}}}\\ {{C_n} ={-} {a_{0,n}}} \end{array}} \right.. $$

In this way, the proposed CMGDM can be expressed as:

(20)$${\varepsilon _r}(\omega )= {\varepsilon _\infty } + \mathop \sum \nolimits_{n = 1}^N (\frac{{{a_{0,n}}}}{{{b_{0,n}} + j\omega }} + \frac{{{C_n}}}{{j\omega }}) + \mathop \sum \nolimits_{m = N + 1}^M \frac{{{a_{0,m}} + j\omega {a_{1,m}}}}{{{b_{0,m}} + j\omega {b_{1,m}} - {\omega ^2}}}, $$

where:

(21)$${C_n} = \left\{ {\begin{array}{{c}} { - {a_{0,n}}\; \; \; \; ,Drude}\\ {0\; \; \; \; \; \; \; \; \; ,others} \end{array}} \right.. $$

This extra term looks artificial, but in fact, it can be incorporated into the ADE quite conveniently. By substituting Eq. (19) into Eq. (2), we have:

(22)$$\sigma \vec{E} + {\varepsilon _0}j\omega \vec{E}\; \left[ {{\varepsilon_\infty } + \left( {\mathop \sum \nolimits_{n = 1}^N \frac{{{a_{0,n}}}}{{{b_{0,n}} + j\omega }} + \mathop \sum \nolimits_{n = 1}^N \frac{{{C_n}}}{{j\omega }}} \right) + \mathop \sum \nolimits_{m = N + 1}^M \frac{{{a_{0,m}} + j\omega {a_{1,m}}}}{{{b_{0,m}} + j\omega {b_{1,m}} - {\omega^2}}}} \right] = \nabla \times \vec{H}. $$

The extra term $\frac{{{C_n}}}{{j\omega }}$ can be transformed into the expression for conductivity as follows:

(23)$$\left( {\sigma + \mathop \sum \nolimits_{n = 1}^N {\varepsilon_0}{C_n}} \right)\vec{E} + {\varepsilon _0}j\omega \vec{E}\left( {{\varepsilon_\infty } + \mathop \sum \nolimits_{n = 1}^N \frac{{{a_{0,n}}}}{{{b_{0,n}} + j\omega }} + \mathop \sum \nolimits_{m = N + 1}^M \frac{{{a_{0,m}} + j\omega {a_{1,m}}}}{{{b_{0,m}} + j\omega {b_{1,m}} - {\omega^2}}}} \right) = \nabla \times \vec{H}. $$

The first term on the left-hand side indicates that the introduced term ${C_n}$ can be integrated into Maxwell’s equation as a polarization conductivity.

3.2 Parameters of the conductive mixed-order generalized dispersive model

The proposed CMGDM does not change the parameters of the 2^nd order dispersion terms (Lorentz, Sellmeier, and Critical Points). But it represents the Debye and Drude model by the 1^st order Padé polynomial and a conductive term. Based on the above derivation, the parameters of the 1^st order CMGDM can be summarized as follows:

Similar to the GDM, the proposed CMGDM results from the Padé approximation of an argument $j\omega $ with real coefficients. Therefore, it does not naturally conform to causal physics, namely the Kramers-Kronig (KK) relations. However, in most applications, the parameters of the GDM can be selected to obey causality. For example, in this paper, we used the D2CP model, which contains a Drude term and two critical point terms, that automatically satisfy KK-consistency.

By reviewing prior D2CP studies, we selected models for Au, Ag, and Al that display wideband accurate fitting. Following the relationship shown in Table 1 (2^nd order) and Table 3 (1^st order), also considering the dispersive models of Eq. (7) and (8) in wavelength form, the parameters of the proposed CMGDM for several noble metals are provided in Table 4.

Table 3. The 1^st order CMGDM parameters for 1^st order dispersion terms

View Table | View all tables in this article

Table 4. Recommended CMGDM parameters for noble metals

View Table | View all tables in this article

As demonstrated in Figs. 1, 2, and 3, the proposed CMGDM yields a very good fit with experimental data over a wide bandwidth. Generally, the refractive index has a relative error below 30% at all wavelengths of interest, which is consistent with other wideband fitting models reported in the literature. Note that through the derivations presented in Section 2.1, the GDM is mathematically equivalent to the D2CP model. Moreover, based on the derivations in Section 3.1, the CMGDM is mathematically equivalent to the traditional GDM. Thus, the proposed CMGDM model has identical fitting accuracy and generality as the traditional GDM method. In the next section, we will discuss the advantage of the CMGDM in terms of fewer unknowns when implemented into an ADE technique.

Fig. 1. The (a) permittivity and (b) refractive index of Au resulting from the proposed CMGDM and experimental measurement [19]. The relative fitting error for the (c) permittivity and (d) refractive index.

Download Full Size | PDF

Fig. 2. The (a) permittivity and (b) refractive index of Ag resulting from the proposed CMGDM and experimental measurement [20]. The relative fitting error for the (c) permittivity and (d) refractive index.

Download Full Size | PDF

Fig. 3. (a) Permittivity and (b) refractive index of Al resulting from the proposed CMGDM and experimental measurement [21]. The relative fitting error for the (c) permittivity and (d) refractive index.

Download Full Size | PDF

3.3 ADE implement of CMGDM model

To implement the D2CP CMGDM model into the DGTD time domain solver, first we rewrite Eq. (2) by taking Eq. (23) into consideration:

(24)$$\left( {\sigma + \mathop \sum \nolimits_{n = 1}^N {\varepsilon_0}{C_n}} \right)\vec{E} + {\varepsilon _0}j\omega \vec{E}\left( {{\varepsilon_\infty } + \mathop \sum \nolimits_{n = 1}^N \frac{{{a_{0,n}}}}{{{b_{0,n}} + j\omega }} + \mathop \sum \nolimits_{m = N + 1}^M \frac{{{a_{0,m}} + j\omega {a_{1,m}}}}{{{b_{0,m}} + j\omega {b_{1,m}} - {\omega^2}}}} \right) = \nabla \times \vec{H}, $$

(25)$$\nabla \times \vec{H} = \left( {\mathop \sum \nolimits_{n = 1}^N {\varepsilon_0}{C_n}} \right)\vec{E} + {\varepsilon _0}{\varepsilon _\infty }j\omega \vec{E} + {\varepsilon _0}j\omega \mathop \sum \nolimits_{n = 1}^N {\vec{P}_n} + {\varepsilon _0}j\omega \mathop \sum \nolimits_{m = N + 1}^M {\vec{P}_m}$$

For D2CP metals, we have $\sigma = 0,\; \; N = 1,\; M = 3$. The 2^nd order polarization vectors ${\vec{P}_{N + 1}}$,…, ${\vec{P}_M}$ are defined the same as in Eq. (12), while the 1^st order polarization ${\vec{P}_n}$ is:

(26)$${\vec{P}_n} = \frac{{{a_{0,n}}}}{{{b_{0,n}} + j\omega }}\vec{E}. $$

In contrast with Eq. (12), however, it possesses only a 1^st order term in $\omega $. Therefore, it does not require an additional auxiliary vector ${\vec{Q}_1}(t )$ as in (13).

By applying an inverse Fourier transform, and considering Eq. (1), (12), (13), (25) and (26), we arrive at:

(27)$$\nabla \times \vec{E} ={-} {\mu _0}{\mu _r}\frac{{\partial \vec{H}(t )}}{{\partial t}}, $$

(28)$$\nabla \times \vec{H} = {\varepsilon _0}C\vec{E} + {\varepsilon _0}{\varepsilon _\infty }\frac{{\partial \vec{E}(t )}}{{\partial t}} + {\varepsilon _0}\frac{{\partial {{\vec{P}}_1}(t )}}{{\partial t}} + {\varepsilon _0}[{{{\vec{Q}}_2}(t )+ {{\vec{Q}}_3}(t )} ]\;,$$

(29)$$\frac{{\partial {{\vec{P}}_1}(t )}}{{\partial t}} = {a_{0,1}}\vec{E}(t )- {b_{0,n}}{\vec{P}_1}(t ), $$

(30)$$\frac{{\partial {{\vec{Q}}_m}(\textrm{t} )}}{{\partial t}} + {b_{1,m}}{\vec{Q}_m}(\textrm{t} )+ {b_{0,m}}{\vec{P}_m}(\textrm{t} )= {a_{0,m}}\vec{E}(t )+ {a_{1,m}}\frac{{\partial \vec{E}(t )}}{{\partial t}},m = 2,3, $$

(31)$${\vec{Q}_m}(t )= \frac{{\partial \vec{P}(t )}}{{\partial t}},\; m = 2,3. $$

The above expressions (27)-(31) can be represented as 7 equations and 7 unknown vectors in order to solve for the quantities ($\vec{E}$, $\vec{H}$, ${\vec{P}_1}$, ${\vec{P}_2}$, ${\vec{P}_3}$, ${\vec{Q}_2}$, ${\vec{Q}_3}$). These are solved in our DGTD formulation by using a fourth-order, four-stage explicit Runge-Kutta method (ERK). The detailed process can be found in Ref. [13].

Note that there are only 7 unknown vectors instead of the 8 derived in Section 2. Therefore, comparing to the conventional GDM method, the proposed CMGDM method has 1/8 less unknowns. Therefore, it is expected to exhibit about a 12.5% higher efficiency and reduced memory consumption. Moreover, as the CMGDM is mathematically equivalent to the GDM method, it maintains the same accuracy.

4. Numerical examples

In this section, we validate the accuracy and efficiency of the proposed algorithm. To accomplish this goal of validating the improvements of the proposed CMGDM model over the conventional GDM, we analyzed the reflection and transmission of a gold slab upon normal incidence. The analytical result can be derived from the well-known closed-formed Fresnel equations. By comparing the numerical results with the analytical solution, we can determine the accuracy and convergence rate of the proposed method and conventional GDM.

Figure 4(a) shows the geometry under investigation. The gold slab (orange) is assumed to be infinite in extent and truncated by periodic boundaries. It has a thickness of 36 nm, sandwiched by air (blue), and is illuminated by a sinusoidally modulated Gaussian pulse with a magnetic current excitation of amplitude ${A_0}(t )$:

(32)$${A_0}(t )={-} \textrm{exp} \left[ { - \frac{{4\mathrm{\pi }{{({t - {t_0}} )}^2}}}{{{t_m}}}} \right]\textrm{cos}[{2\mathrm{\pi }{f_m}t} ], $$

where ${t_0}$ denotes the reference temporal points, ${t_m}$ is a spectral bandwidth related parameter, and ${f_m}$ is the frequency of the modulated signal. In the direction of illumination, the computational domain is truncated by two wave ports. The transmission and reflection of the gold-air boundaries can be theoretically calculated by Fresnel equations. In order to better highlight the efficiency improvement of the proposed CMGDM method, the studied gold slab is assumed to be thick, occupying 90% of the computational domain’s total volume.

Fig. 4. (a) Computational domain of the numerical example, which is truncated by periodic boundaries and a pair of wave ports. (b)Transient response of a gold slab under x-polarized normally incident illumination. (c) Transmission and reflection performance in wavelengths. (d) Convergence plot showing accuracy improvement along with the refined mesh size.

Download Full Size | PDF

Figure 4(b) shows the transient response of the gold slab when subject to an x-polarized normally incident plane wave. The reflection and transmission performance of the gold slab can be extracted through Fourier transformation, as shown in Fig. 4(c). As can be seen, the results of the proposed CMGDM and the conventional GDM match very well. Figure 4(d) shows convergence plots of CMGDM and GDM that have been generated by comparing them with the analytical result. This provides a means for quantifying the degree of accuracy improvement against refinements in the mesh [4,14–18]. The convergence plots for CMGDM and conventional GDM are almost identical as they are mathematically equivalent with each other.

For the case where $\lambda /h = 10$ where h is the maximum value of the mesh length, the GDM method requires 32.0 MB of memory and 15.1 seconds of computation time, while the proposed CMGDM method requires only 29.1 MB memory and 13.5 seconds. This means that the proposed CMGDM reduces memory consumption by 9% while simultaneously accelerating the computation by 10.6% compared to the conventional GDM method. In order to reduce the environmental perturbation and artificial error, we increased the computational complexity by using an over-refined mesh with $\lambda /h = 100$. For the refined mesh, the CMGDM outperforms the conventional GDM with memory consumption reduced by 9% from 879 MB to 800 MB. Moreover, the computational time is accelerated by 11% from 8542 seconds to 7691 seconds. The efficiency improvement comes from the lower differential order of CMGDM, which requires fewer unknowns. Moreover, it should be pointed out that the efficiency improvement of CMGDM is achieved without sacrificing accuracy.

All numerical simulations were performed on a laptop with Dual Intel E5-2630 CPU of 2.4 GHz, 32 threads, and 144 GB of memory. The algorithm has been fully parallelized with all 32 threads using OpenMP. The comparison was made with the DGTD GDM simulation tool that has been released to public [22].

5. Conclusions

A CMGDM model was proposed for accelerating the simulation of dispersive metals in the optical regime. It is composed of a series of mixed-order Padé polynomials with an extra conductive term. The parameters of the GDM and CMGDM were derived from selected D2CP fitting datasets. When compared to the conventional GDM model, the CMGDM has lower differential order, and thus fewer auxiliary unknown vectors. Consequently, the CMGDM method provides higher computational efficiency, lower memory consumption, while maintaining the same accuracy.

Funding

Pennsylvania State University (MRSEC); National Science Foundation (DMR-1420620); Defense Advanced Research Projects Agency (HR00111720032).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

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

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

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

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

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

6. P. G. Etchegoin, E. C. Le Ru, and M. Meyer, “An analytic model for the optical properties of gold,” The Journal of Chemical Physics 125(16), 164705 (2006). [CrossRef]

7. A. Vial, T. Laroche, M. Dridi, and L. Le Cunff, “A new model of dispersion for metals leading to a more accurate modeling of plasmonic structures using the FDTD method,” Appl. Phys. A 103(3), 849–853 (2011). [CrossRef]

8. M. Garcia-Vergara, G. Demésy, and F. Zolla, “Extracting an accurate model for permittivity from experimental data: hunting complex poles from the real line,” Opt. Lett. 42(6), 1145 (2017). [CrossRef]

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

10. L. J. Prokopeva, J. Trieschmann, T. A. Klar, and A. V. Kildishev, “Numerical modeling of active plasmonic metamaterials,” in Optical Complex Systems: OCS11 (International Society for Optics and Photonics, 2011), 8172, p. 81720B.

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

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

13. W. Mai, S. D. Campbell, E. B. Whiting, L. Kang, P. L. Werner, Y. Chen, and D. H. Werner, “Prismatic discontinuous Galerkin time domain method with an integrated generalized dispersion model for efficient optical metasurface analysis,” Opt. Mater. Express 10(10), 2542–2559 (2020). [CrossRef]

14. D. Barchiesi and T. Grosges, “Fitting the optical constants of gold, silver, chromium, titanium, and aluminum in the visible bandwidth,” J. Nanophotonics 8(1), 083097 (2014). [CrossRef]

15. 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).

16. A. D. Rakić, A. B. Djurišić, J. M. Elazar, and M. L. Majewski, “Optical properties of metallic films for vertical-cavity optoelectronic devices,” Appl. Opt. 37(22), 5271–5283 (1998). [CrossRef]

17. H. S. Sehmi, W. Langbein, and E. A. Muljarov, “Optimizing the Drude-Lorentz model for material permittivity: Method, program, and examples for gold, silver, and copper,” Phys. Rev. B 95(11), 115444 (2017). [CrossRef]

18. R. F. Ando, A. Tuniz, J. Kobelke, and M. A. Schmidt, “Analysis of nanogap-induced spectral blue-shifts of plasmons on fiber-integrated gold, silver and copper nanowires,” Opt. Mater. Express 7(5), 1486–1495 (2017). [CrossRef]

19. P. B. Johnson and R. W. Christy, “Optical Constants of the Noble Metals,” Phys. Rev. B 6(12), 4370–4379 (1972). [CrossRef]

20. P. Winsemius, F. F. van Kampen, H. P. Lengkeek, and C. G. van Went, “Temperature dependence of the optical properties of Au, Ag and Cu,” J. Phys. F: Met. Phys. 6(8), 1583–1606 (1976). [CrossRef]

21. A. D. Rakić, “Algorithm for the determination of intrinsic optical constants of metal films: application to aluminum,” Appl. Opt. 34(22), 4755 (1995). [CrossRef]

22. Wending Mai and Douglas H. Werner, Prism-DGTD with GDM to Analyze Pixelized Metasurfaces (Zenodo, 2020).

		$a_{0}$	$a_{1}$	$b_{0}$	$b_{1}$
Au	1	1.7351e32	0	0	1.2991e14
ɛ_∞= 1.54	2	3.5985e31	7.1982e15	1.7045e31	1.9828e15
	3	6.8279e31	9.0162e15	3.6750e31	3.5541e15
Ag	1	1.7636e32	0	0	9.1269e13
ɛ_∞= 1.4447	2	2.2018e32	3.377e14	9.5275e31	1.0354e16
	3	2.9375e30	3.1258e15	3.8730e31	1.0824e15
Al	1	4.2428e32	0	0	2.2876e14
ɛ_∞= 1.0	2	5.0791e31	1.1631e16	5.2582e30	6.5734e14
	3	3.9424e31	1.0722e16	9.2290e30	3.5462e15

	Debye	Drude
$a_{0}$	$C_{D e} / τ_{D e}$	$- ω_{D r}^{2} / τ_{D r}$
$b_{0}$	$1 / τ_{D e}$	$τ_{D r}$
C	0	$ω_{D r}^{2} / τ_{D r}$

		$a_{0}$	$a_{1}$	$b_{0}$	$b_{1}$	C
Au ɛ_∞= 1.54	1	−1.3356e18	/	1.2991e14	/	1.3356e18
	2	3.5985e31	7.1982e15	1.7045e31	1.9828e15	/
	3	6.8279e31	9.0162e15	3.6750e31	3.5541e15	/
Ag ɛ_∞= 1.4447	1	−1.9323e18	/	9.1269e13	/	1.9323e18
	2	2.2018e32	3.377e14	9.5275e31	1.0354e16	/
	3	2.9375e30	3.1258e15	3.8730e31	1.0824e15	/
Al ɛ_∞= 1.0	1	−1.8547e18	/	2.2876e14	/	1.8547e18
	2	5.0791e31	1.1631e16	5.2582e30	6.5734e14	/
	3	3.9424e31	1.0722e16	9.2290e30	3.5462e15	/

		$a_{0}$	$a_{1}$	$b_{0}$	$b_{1}$
Au	1	1.7351e32	0	0	1.2991e14
ɛ_∞= 1.54	2	3.5985e31	7.1982e15	1.7045e31	1.9828e15
	3	6.8279e31	9.0162e15	3.6750e31	3.5541e15
Ag	1	1.7636e32	0	0	9.1269e13
ɛ_∞= 1.4447	2	2.2018e32	3.377e14	9.5275e31	1.0354e16
	3	2.9375e30	3.1258e15	3.8730e31	1.0824e15
Al	1	4.2428e32	0	0	2.2876e14
ɛ_∞= 1.0	2	5.0791e31	1.1631e16	5.2582e30	6.5734e14
	3	3.9424e31	1.0722e16	9.2290e30	3.5462e15

	Debye	Drude
$a_{0}$	$C_{D e} / τ_{D e}$	$- ω_{D r}^{2} / τ_{D r}$
$b_{0}$	$1 / τ_{D e}$	$τ_{D r}$
C	0	$ω_{D r}^{2} / τ_{D r}$

Conductive mixed-order generalized dispersion model for noble metals in the optical regime

Abstract

1. Introduction

2. Conventional D2CP GDM model

2.1 Parameters of GDM model

2.2 ADE implement of conventional GDM model

3. Conductive mixed-order GDM model

3.1 Derivation

3.2 Parameters of the conductive mixed-order generalized dispersive model

3.3 ADE implement of CMGDM model

4. Numerical examples

5. Conclusions

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (4)

Tables (4)

Equations (32)

Optics Express