The discontinuous Galerkin time domain (DGTD) method and its recent flavor, the continuous-discontinuous Galerkin time domain (CDGTD) method, have been extensively applied to simulations in the radio frequency (RF) and microwave (MW) regimes due to their inherent ability to efficiently model multiscale problems. We propose to extend the CDGTD method to nanophotonics while exploiting its advantages which have already been established in the RF and MW regimes, such as domain decomposition, non-conformal meshing, high-order elements, and hp-refinement. However, at optical frequencies many materials are highly dispersive, so the modeling of nanophotonic devices requires accurate handling of different dielectric functions, including those of plasmonic elements, dielectrics, and tunable materials. In this paper, we propose a CDGTD method that incorporates a generalized dispersive material (GDM) model which is an efficient way to implement a wide range of optical dispersion models with a universal analytic function. Physics-based dispersion models, such as the Drude, Debye, Lorentz, and critical points as well as more complicated behavior founded on ab-initio principles can all be obtained as special cases of the universal GDM approach. The accuracy and convergence of this GDM-incorporated CDGTD are verified by numerical examples. The CDGTD method, equipped with the GDM model, paves the way to the efficient design and optimization of large scale photonic devices with a diverse range of optical dispersive materials.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
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  models along with some classical empirical models (e.g., the Cauchy and Cole-Cole models) . 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  was successfully implemented in DGTD simulations of graphene , 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 ADE1,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 a0,1, b0,1 (3) that can be obtained from Table 1.
2.2 CDGTD method incorporating the GDM model
The first-order Maxwell’s equations with the polarization currents can be expressed as follows21–27]. The polarization currents, Jm 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 kth subdomain to be written as22,23].
Next, by applying central differences to (6) and (8) at and to (7) at we can obtain the system of linear equations as follows,
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 . 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 and being the conversion formulas to their original values in rad/s, respectively.
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  and listed in Table 2. The incident field is generated by a modulated Gaussian pulse
For the detailed implementation of the plane wave excitation for the vector basis functions, please refer to . 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 , 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 asFig. 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).
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 . 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  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 and fields, instead of and 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  as shown in Fig. 3(b). The outer region and the perfectly matched layer (PML) are meshed with high-order brick elements (not shown).
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.
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.
DARPA/DSO Extreme Optics and Imaging (EXTREME) Program, Award # HR00111720032.
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]
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).
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]