We present a modified formulation of the Finite-Difference Time-Domain (FDTD) technique that facilitates the accurate modeling of curved plasmonic interfaces. These interfaces appear in structures of interest for the design of optical metamaterials, such as arrays of plasmonic nanorods. Our approach uses the standard rectangular FDTD mesh and tensor effective permittivities for the interface cells, implicitly enforcing field boundary conditions, and is readily applicable to thin curved dispersive layers. We demonstrate the accuracy and effectiveness of our approach with the periodic analysis of a silver nanorod array and the computation of scattering parameters from a thin dispersive ring in a waveguide.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Various geometries of interest involve curved interfaces between dielectrics and metals. These include periodic split ring resonators (SRR) and plasmonic nanorods, which have been employed in a wide range of applications related to optical metamaterials [1–5]. Another important geometry is the plasmonic ring resonator which has been recently used to control the transmission of surface plasmon polaritons (SPPs) [6, 7]. These motivate the analysis of the wave propagation properties of these structures with full-wave numerical techniques such as the Finite-Difference Time-Domain (FDTD) method. The FDTD method is a popular numerical tool for photonics on account of its ability for handling complex geometries and dispersive media, and its broadband characteristics. However, the method suffers from a serious accuracy degradation due to pronounced staircasing errors when approximating curved interfaces by its conventional rectangular grid [8, 9]. Notably, the staircasing errors become more severe in plasmonic structures, as surface plasmons are strongly localized at metal-dielectric interfaces, where staircasing is applied.
Previous research on these problems is mainly based on two different approaches: modifying effective permittivities to account for the interface cells [10–13] and unstructured mesh techniques [14, 15]. A scalar effective permittivity technique was adopted in [10, 11] and extended to dispersive media by a z-transform and an auxiliary differential equation (ADE) method, respectively. However, the use of scalar effective parameters actually violates the field boundary conditions at the interface, and the resultant fourth-order time derivatives in  need significant memory and execution time. Another effective permittivity theory of interest is the subpixel smoothing technique , which employed a tensor effective permittivity. This technique was extended to sloped dispersive interfaces in  by splitting the electric fields into three auxiliary variables, instead of one, which greatly increases memory and computational requirements. By employing a coordinate rotation of the effective permittivity tensor,  reduced the three variables back to two. Moreover,  introduced a triangular mesh FDTD scheme in 2-D, which is provably stable and rapidly convergent. However, the time step limit is severely constrained by the smallest element of the triangular mesh. Along similar lines,  derived a linearized non-local dispersion model using the Discontinuous Galerkin Time-domain (DGTD) method. The advantage of these techniques is that the unstructured triangular mesh naturally achieves better conformity to curved interface structures than the rectangular FDTD mesh.
In this paper, we present a different approach that strikes a good balance between accuracy and efficiency, by keeping the structured mesh of FDTD and implicitly enforcing field boundary conditions at the interface Yee cells. This is achieved by formulating a tensor ADE-FDTD technique, extending the tensor FDTD method of [17, 18] from simple lossy dielectrics to dispersive media of Lorentz type dispersion. Through further derivation, the proposed method can also handle thin subcellular curved sheets in a Yee cell. Both [19, 20] independently introduced effective permittivities for dielectric thin sheet inclusions by enforcing boundary conditions at material interfaces. However, when using a structured FDTD mesh, these methods only work with planar interfaces. In this paper, we introduce an effective tensor permittivity to account for the boundary conditions at the two interfaces of thin sheets, thus making it possible to embed arbitrarily inclined and curved thin sheets in a Yee cell.
The rest of the paper is structured as follows: in section II, the effective tensor permittivity for sloped interfaces and thin sheets is derived. In section III, the tensor ADE-FDTD update scheme for Lorentz media is derived. In section IV, periodic silver nanorods are simulated to demonstrate the proposed method’s effectiveness at handling sloped dispersive interfaces. As well, a PEC cavity with dielectric thin sheets is modeled to validate our thin sheet model. Finally, a thin ring of a Lorentz medium is simulated to show the efficiency of the tensor ADE-FDTD for dispersive, curved thin sheets. Section V summarizes our contributions.
2. Effective tensor permittivity for sloped interfaces and thin sheets
The adopted tensor effective permittivity model was proposed in  for handling simple, dielectric interfaces. Then the technique was improved by , which adjusted the tensor effective permittivity to field polarization and was extended to simple, lossy media interfaces in , via a z-transform method. In this paper, we adopt the tensor model of  and consider the 2-D case of Fig. 1(a), where fields (Ex, Ey, Hz) are sought for. The FDTD update scheme will be derived based on the integral form of Maxwell’s equations. The FDTD update equation for , for instance, can be derived by formulating a loop (Fig. 1(b)), which has length of along the direction, as in the direction both fields and modeled geometries are considered to be invariant. Thus, the flux density can be updated as the following expression:
Note that and are at the same grid point. If the medium inside the face is homogeneous, the determination of from only requires the medium permittivity. However, if the medium is inhomogeneous, boundary conditions must be enforced. In the following, we will first derive the effective tensor permittivity for sloped interfaces and then extend it to sloped thin sheets.
Considering two Yee cells centered at the grid point , the background dielectric medium ε1 is interfaced with a dispersive medium , as in Fig. 2(a). In the following, we discuss how to find the tensor effective permittivity for the update of the electric field node .
In the middle face centered at and between the two magnetic nodes, we consider the electric field in the background medium and the electric field in the dispersive medium. The relationship between and depends on the interface normal and must satisfy the following boundary conditions:Eq. (2) and Eq. (3) we obtain: Fig. 2(b)). We consider the electricfield above the thin sheet, ; the electric field inside the sheet, and the electric field below the sheet, . By applying boundary conditions of Eqs. (2) and (3) to pairs , and , we obtain:
The matrix is the same for the two pairs of fields as the interface normal is the same and introduces no difference to the matrix in Eq. (4). As a result, we can derive that . By replacing with , the electric field relation of Eq. (4) can be directly adopted in the thin sheet model. Notably, sloped interfaces and thin sheets share the same relationship between the electric field in the dispersive medium and that in the background medium. Once we obtain the relationship between the electric fields and , we can derive the effective tensor permittivity. The electric flux through a face of a Yee cell can be expressed in terms of either D or E as:
We can apply the above equation to to a face similar to Fig. 1(b). Figures. 3(a) and 3(b) show the corresponding faces for the sloped interface and the thin sheet. The part of the face in the dispersive medium is depicted as gray, and their widths are indicated as g for both cases, which are also marked in Fig. 2. The flux density can then be expressed as follows:
Hence:Eq. (4) we obtain:
Note that in Eq. (9), the flux density Dy is only dependent on the electric field in medium 1. However, the electric field in the direction also contributes to Dy, which means that the electric field in the direction cannot be directly derived by Dy. Here, we introduce another flux component , at the same grid point as . By repeating the above procedure via Eqs. (6)–(9) for , the relationship between the flux density D at the grid point and the electric field is formulated as:
3. ADE-FDTD update scheme
In the following, we will first consider the medium in region 2 as a simple dielectric medium and give the normal update scheme. Then, we will consider the medium in region 2 as a Lorentz dispersive medium, and finally give the ADE-FDTD update scheme.
3.1. Electric field update
In this section, we focus on the update scheme for the electric field using the derived tensor effective permittivity in Eq. (10). We present the update of the electric field node . Let us first consider the medium in region 2 as a simple dielectric medium. In every FDTD time step, we can readily update D having the magnetic field at time step n+1/2. Then, D can be used to derive from:
To derive from the above equation, we need all components of D at the interface grid node. is found by the normal update equation of Eq. (1) and can be found by interpolating the adjacent Dx components:
Once is derived from the flux density using , we can use Eq. (4) to find . Then, we need to check whether the electric field node is in the dielectric host or in the dispersive medium. If it is within the dielectric, is updated by the y-component of ; otherwise is updated by the y-component of . Note that the y-components of both and need to be stored and will be used in the update of the magnetic field.
Then, we consider the ADE-FDTD update scheme for cases when the region 2 is filled with a dispersive medium. The following Lorentz dispersion model is adopted to describe the permittivity of the dispersive medium:Eq. (11) to directly derive the electric field in the dielectric host. To solve this problem, we introduce an auxiliary vector as:
Then, Eq. (10) is cast in the following form:
Mapping Eq. (14) to the time-domain and discretizing it with second-order finite-differences yields update equations of the form:
The coefficients are given by:
Then we need to derive the electric field in the dispersive medium via Eq. (4). However, cannot be directly derived, due to its frequency-dependent relation to . Here we introduce another auxiliary vector as:
Then, Eq. (4) becomes:
Mapping the auxiliary vector to the time-domain and discretizing it with second-order finite-differences yields update equations of the following form (where the coefficients and are the same as in Eq. (17)):
Finally, we check the location of the electric node to decide whether it is updated by or by . We also need to store both and for the update of the magnetic field.
3.2. Magnetic field update
Once all the electric field components have been computed, the magnetic field updates can be derived by contour path integration along the perimeter of a standard cell, as shown in Fig. 4. The lengths inside the dispersive medium on the left and bottom cell edges are marked as ly and lx for both models. Then, the line integration of the electric field along the left cell edge is expressed as:
3.3. Implementation of tensor ADE-FDTD for a Lorentz dispersive medium
Before beginning the FDTD time loop, the coefficients of Eq. (17) are determined and stored. Then, at every cell on the curved interface or sloped thin sheet, the matrices ,, and are stored for the updates of the electric field.
For electric field updates: 1) Compute Dx and Dy at every interface grid node via the standard FDTD update equations and interpolate, at intervals, to get the Dy and Dx at each other’s grid points. 2) Compute at every interface cell according to Eq. (18), then compute at every interface cell via Eq. (22). 3) Update the electric field by or based on the location of the corresponding electric field nodes. 4) Update the auxiliary vectors P and Q using Eq. (16) and Eq. (21).
4. Numerical results
In this section, we demonstrate the proposed tensor ADE-FDTD for sloped interfaces by modelling a one-dimensional array of plasmonic cylinders. Then, we apply the thin sheet version of the method to model a thin Lorentzian ring in a waveguide.
4.1. 1-D periodic silver nanorods
We will consider an array of silver nanorods (Fig. 5), an example also treated in [10, 14]. The one dimensional periodic structure is modeled by Bloch’s periodic boundary conditions. In our two dimensional TE computational domain, Ex and Hz are involved in the update equations and the periodic boundary conditions, as tangential to the direction of periodicity (here defined by the y-axis):Eq. (13). A Gaussian point source is used to excite the Bloch modes of the periodic structure. A spectral analysis of the fields provides the corresponding resonant frequencies. Fig. 6 presents the Fourier transform of the electric field at a sampling point within the unit cell for , for the tensor ADE-FDTD and a standard, staircased FDTD with different cells of area . For the staircased FDTD, we assume that the electric field nodes are filled with the dispersive medium if the nodes are located inside the cylinder. Otherwise the electric field nodes are updated by the permittivity of the dielectric host (vacuum here). The results of the staircased FDTD change significantly as the mesh is refined. Notably, staircasing contributes to spurious modes as well, which were also observed in . On the other hand, the tensor ADE-FDTD method leads to a correct resonant frequency for cell size nm; standard FDTD needs nm or smaller cell size to achieve a similar accuracy.
The frequencies of the first five modes are shown in Table 1, along with those found by applying the isotropic effective permittivities (EPs) method of . The results of the proposed tensor ADE-FDTD method converge well as the mesh is refined, unlike the results of the staircased FDTD which exhibit oscillations with mesh refinement. Then, we compare the resonant patterns of the first five modes using both methods. The accurately resolved modal patterns of the magnetic field along with the electric field vectors are plotted in Fig. 7, using the tensor ADE-FDTD method. Similar graphs are generated using the staircased FDTD method in Fig. 8. The cell size is chosen as 2.5 nm for both methods. The patterns found with the proposed method are clear and symmetrical. In comparison, the staircased FDTD fails to obtain clear modal field patterns, especially for the last three modes, although their corresponding resonant frequencies are close to those of the proposed tensor ADE-FDTD method. Then, we check the convergence rates of the two methods. We compute the root-mean-square (RMS) error of the frequencies of the first two modes, with the following equation:Fig. 9. Notably, the proposed method shows a faster and monotonic convergence. However, the staircased FDTD method converges slower and the results oscillate as the mesh is refined. The execution time and RMS error for the first two modes are compared in Table 2. Finally, the Bloch diagram is shown and compared with that of  in Fig. 10. For the first two modes, the two methods agree well. However, the results of  deviate from the results ofthe tensor ADE-FDTD method when higher modes are considered. This is likely due to the implicit enforcement of field boundary conditions at the metal-dielectric interface in the proposed method, in contrast to the scalar effective permittivity based method of .
4.2. Sloped thin dispersive sheet
In this section, we present two case studies, where the proposed tensor ADE-FDTD technique is applied to sloped thin sheets. Our method is first evaluated in a cavity geometry and then applied to compute the scattering parameters of a thin dispersive ring in a waveguide.
4.2.1. Verification of thin sheet model with simple dielectric media
The accuracy of our proposed method was verified by simulating a 50 cm × 2.5 cm cavity with PEC walls loaded with parallel, sloped thin sheets. The layout of the simulation is shown in Fig. 11(a). The dielectric thin sheets (, μ0) are located in free space (ϵ0, μ0). Note that in this case, the permittivity in region 2 is a constant, so the effective tensor permittivity in Eq. (10) is a simple matrix. Hence, the normal update scheme via Eq. (11) can be adopted. The cavity was excited by a Gaussian point source with a bandwidth of 3GHz. To verify the accuracy of the proposed method, we compared the resonant frequency of the structure computed by the tensor FDTD and standard FDTD with a very fine mesh (reference FDTD). In the reference FDTD simulation, the structure was discretized by cells much smaller than the sheet width, with a cell size mm; the same results were obtained with even smaller cells. The cell size of tensor FDTD was chosen as mm, where vp is the speed of light inside the dielectric sheets. The width of the thin sheet is mm and the inclination angle is . Fig. 11(b) illustrates the comparison between tensor FDTD and reference FDTD with respect to the spectrum ofthe electric field within the cavity. The results agree well. We also test the accuracy of the method with varying inclination angles of the thin sheets. The sheets with a fixed width of w = 1 mm are tilted from 0° to 45° with a step of . The errors (compared to the reference FDTD) of three modes (TE, TE and TE) are shown in Fig. 12(a). The accuracy remains high as θ varies, which proves the effectiveness of our method for arbitrary angles. Finally, we checked the error variation with sheet width. We fixed the sheet angle to and changed the sheet width from 0.5 mm to 3 mm with a step size of 0.5 mm. Fig. 12(b) illustrates how the error of the three modes varies with sheet width. The results show that the error is very small when the sheet is thin. However, the accuracy decreases as the sheet width increases. This phenomenon is due to the assumption that the electric field in region 2 is constant. The variation of the electric field inside the thin sheet can be ignored when the sheet is sufficiently thin. Then, Eq. (5) is satisfied. As the sheet becomes thicker, the electric field differs across the two interfaces, thus Eq. (5) cannot be satisfied. As a result, the accuracy of the method degrades in this case.
4.2.2. Thin Lorentz medium ring in a waveguide
Finally, let us consider a thin ring filled with a Lorentz medium, inside a waveguide (Fig. 13). The ring has an inner radius of nm and an outer radius of nm. For the Lorentz medium parameters in Eq. (13), we use rad/s, rad/s and Hz. The ring is placed in a parallel plate waveguide, excited by a plane wave with maximum frequency of Hz. The scattering parameters of the ring are computed by probing the vertical electric field in front of and behind the ring. With the tensor ADE-FDTD method, we introduce cells much larger than the width of the ring, which are marked with red squares in Fig.13. The S-parameters are computed and compared to the standard, staircased FDTD simulation in fine mesh (Fig. 14). The tensor ADE-FDTD method adopts a cell mm, where vp is the wave velocity in vacuum. The standard FDTD converges well when the cell size is smaller than 0.0222 nm. A good agreement between the proposed method and the reference FDTD is observed throughout the simulated frequency range. Finally, the execution times are compared in Table 3. Our method requires less than 10 seconds to compute the scattering parameters, where as the staircased FDTD requires two orders of magnitude greater execution time for the same level of accuracy.
A generalized tensor ADE-FDTD method for sloped dispersive interfaces and thin sheets has been presented. The improvement in accuracy for sloped dispersive interfaces is validated by computing the plasmonic modes of a one-dimensional silver nanorods array. Compared with the conventional staircased FDTD, our proposed method shows faster and more stable convergence; plasmonic mode patterns are also well resolved. For the sloped thin dispersive sheet, we have compared the proposed method with an ultra-fine mesh standard FDTD simulation, for the scattering parameters of a thin dispersive ring loaded waveguide. We have shown that our method uses much less time to extract accurate results. The proposed method shows great versatility for arbitrarily shaped interfaces and thin sheets. Hence, it is a useful technique for computational electromagnetics and optics.
Natural Sciences and Engineering Research Council of Canada (NSERC) (453936).
1. A. Alú and N. Engheta, “Optical nanotransmission lines: synthesis of planar left-handed metamaterials in the infrared and visible regimes,” J. Opt. Soc. Am. B 23, 573–583 (2006). [CrossRef]
2. N. Engheta, A. Salandrino, and A. Alú, “Circuit Elements at Optical Frequencies: Nanoinductors, Nanocapacitors, and Nanoresistors,” Phys. Rev. Lett. 95, 095504 (2011). [CrossRef]
4. N. Wongkasem, A. Akyurtlu, K. A. Marx, Q. Dong, J. Li, and W. D. Goodhue, “Development of Chiral Negative Refractive Index Metamaterials for the Terahertz Frequency Regime,” IEEE Trans. Antennas Propag. 55, 3052–3062 (2007). [CrossRef]
5. A. K. Azad, A. J. Taylor, E. Smirnova, and J. F. O’Hara, “Characterization and analysis of terahertz metamaterials based on rectangular split-ring resonators,” Appl. Phys. Lett. 92, 011119 (2008). [CrossRef]
6. T.-B. Wang, X.-W. Wen, C.-P. Yin, and H.-Z. Wang, “The transmission characteristics of surface plasmon polaritons in ring resonator,” Opt. Express 17, 24096 (2009). [CrossRef]
8. A. Cangellaris and D. Wright, “Analysis of the numerical error caused by the stair-stepped approximation of a conducting boundary in FDTD simulations of electromagnetic phenomena,” IEEE Trans. Antennas Propag. 39, 1518–1525 (1991). [CrossRef]
9. S. Dey and R. Mittra, “A locally conformal finite-difference time-domain (FDTD) algorithm for modeling three-dimensional perfectly conducting objects,” IEEE Microw. Guided Wave Lett. 7, 273–275 (1997). [CrossRef]
10. Y. Zhao and Y. Hao, “Finite-Difference Time-Domain Study of Guided Modes in Nano-Plasmonic Waveguides,” IEEE Trans. Antennas Propag. 55, 3070–3077 (2007). [CrossRef]
13. J. Liu, M. Brio, and J. V. Moloney, “Subpixel smoothing finite-difference time-domain method for material interface between dielectric and dispersive media,” Opt. Lett. 37, 4802–4804 (2012). [CrossRef] [PubMed]
14. Y. Liu, C. D. Sarris, and G. V. Eleftheriades, “Triangular-Mesh-Based FDTD Analysis of Two-Dimensional Plasmonic Structures Supporting Backward Waves at Optical Frequencies,” J. Light. Technol. 25, 938–945 (2007). [CrossRef]
15. N. Schmitt, C. Scheid, S. Lanteri, A. Moreau, and J. Viquerat, “A DGTD method for the numerical modeling of the interaction of light with nanometer scale metallic structures taking into account non-local dispersion effects,” J. Comput. Phys. 316, 396–415 (2016). [CrossRef]
16. A. Farjadpour, D. Roundy, A. Rodriguez, M. Ibanescu, P. Bermel, J. D. Joannopoulos, S. G. Johnson, and G. W. Burr, “Improving accuracy by subpixel smoothing in the finite-difference time domain,” Opt. Lett. 31, 2972–2974 (2006). [CrossRef] [PubMed]
17. J. Nadobny, D. Sullivan, W. Wlodarczyk, P. Deuflhard, and P. Wust, “A 3-D tensor FDTD formulation for treatment of sloped interfaces in electrically inhomogeneous media,” IEEE Trans. Antennas Propag. 51, 1760–1770 (2003). [CrossRef]
18. J. Nadobny, D. Sullivan, and P. Wust, “A General Three-Dimensional Tensor FDTD-Formulation for Electrically Inhomogeneous Lossy Media Using the Z-Transform,” IEEE Trans. Antennas Propag. 56, 1027–1040 (2008). [CrossRef]
19. J. Maloney and G. Smith, “The use of surface impedance concepts in the finite-difference time-domain method,” IEEE Trans. Antennas Propag. 40, 38–48 (1992). [CrossRef]
20. K.-P. Hwang and A. Cangellaris, “Effective permittivities for second-order accurate FDTD equations at dielectric interfaces,” IEEE Microw. Wirel. Compon. Lett. 11, 158–160 (2001). [CrossRef]
21. J.-Y. Lee and N.-H. Myung, “Locally tensor conformal FDTD method for modeling arbitrary dielectric surfaces,” Microw. Opt. Technol. Lett. 23, 245–249 (1999). [CrossRef]