We have extended the contour-path effective-permittivity (CP-EP) finite-difference time-domain (FDTD) algorithm by A. Mohammadi et al., Opt. Express 13, 10367 (2005), to linear dispersive materials using the Z-transform formalism. We test our method against staircasing and the exact solution for plasmon spectra of metal nanoparticles. We show that the dispersive contour-path (DCP) approach yields better results than staircasing, especially for the cancellation of spurious resonances.
© 2008 Optical Society of America
The increased popularity of the finite-difference time-domain (FDTD) method as a modeling tool in micro- and nanophotonics  has recently renewed the interest in developing better algorithms to handle material interfaces [2, 3, 4, 5, 6, 7, 8, 9]. In particular, much attention has been devoted to the implementation of effective tensor permittivities to preserve the original Yee meshing scheme [2, 7, 8, 9], which is one of the major advantages of FDTD over other domain discretization methods. Moreover, the advent of left-handed metamaterials [10, 11] and plasmonics [12, 13, 14, 15] has focused the efforts on extending those algorithms to linear dispersive materials [16, 17, 18]. Mohammadi and Agio  have considered only flat metal-dielectric interfaces aligned with the FDTD mesh. Zhao and Hao  have implemented dispersion assuming a diagonal effective permittivity tensor  and the auxiliary differential equation method , obtaining a fourth-order time-stepping scheme that increases both memory and computation time. Deinega and Valuev  have extended the nondiagonal effective permittivity tensor [2, 7, 9] by splitting the electric field in three auxiliary variables such that three, instead of one, auxiliary differential equations are needed. Furthermore, the nondiagonal terms require space interpolation like in the nondispersive case.
In this work, instead of using an effective permittivity tensor, we generalize the contour-path effective-permittivity (CP-EP)  to metal-dielectric interfaces of arbitrary shape. We name this approach dispersive contour-path (DCP) FDTD. By carrying out the implementation of dispersion using the Z-transform  we also show that for flat interfaces the method is equivalent to . Contrary to the effective permittivity tensor, the DCP FDTD algorithm neither requires additional storage nor additional operations in the FDTD algorithm, except for a pre-processing step before time marching the electromagnetic field.
We validate the improvements of our method by computing the plasmon spectra of metal nanoparticles (NPs) and compare the result with the staircase FDTD algorithm  and with exact solutions. We also discuss various issues and limitations of FDTD in modeling metal nanostructures.
2. Dispersive contour-path method
We consider a dispersive nonmagnetic medium with a Drude dielectric function ε 2(ω)=ε ∞-ω 2 p/(ω(ω+iγ)) embedded in a dielectric background ε 1. Here ω p is the plasma frequency and γ is the damping frequency . Following  and  we set both the vacuum permittivity εo and the vacuum permeability µo to 1 and start from Ampère and Faraday laws
and from the constitutive relation using Z-transform
The subscript 2 in the previous equations refers to quantities inside the dispersive material, while the subscript 1 will refer to the same quantities in the dielectric one. We do not use the vector notation because we are considering either the x or the y component on a two-dimensional FDTD mesh. Since flat interfaces parallel or perpendicular to the mesh have been already discussed in the context of surface plasmon polaritons , we focus our attention to the generalization of the DCP approach to partially filled cells with tilted and curved interfaces. Because the treatment of a tilted interface is not fundamentally different from a curved one , without loss of generality we directly consider the situations depicted in Fig. 1.
Figure 1(a) represents a cell where both integration lines cross the interface. Ampère law in Eq. (1) is implemented across the interface by applying the projection technique described in . Using Eq. (2) and defining an average dielectric displacement D in place of D 1 and D 2 we find
where d is the line filling ratio, assuming that the mesh pitch Δ is equal to one, and m is the projection of the interface normal m along the field component. The previous equation can be simplified by grouping the terms containing E 2 and by recognizing that ε 1 and ε ∞ give the effective permittivity ε ‖,m=ε ∞ d+(1-d)[ε ∞ m 2+ε 1(1-m 2)] for nondispersive materials .
We thus write
where c ‖,m=d+(1-d)m 2 derives from the auxiliary term S 2. Notice that it is independent of the dispersion model. The integration in Faraday law is treated in a similar manner by defining an average electric field E
where f is the line filling ratio, assuming that the mesh pitch Δ is equal to one, and n is the projection of the interface normal n along the field component. Using Eq. (2) and grouping the terms containing E 2 we obtain
After taking advantage of Eq. (10) to define an auxiliary term S
Notice that although the coefficients in Eq. (3) depend on the dispersion model, the modifications in Eq. (13) brought by the DCP are general. Furthermore, when ε 2 is nondispersive, the DCP algorithm reduces to the CP-EP FDTD method .
While in the derivation of the CP-EPs the procedure for the FDTD cell in Fig. 1(b) is equal to that for the cell in Fig. 1(a), except that ε ‖,m and ε ⊥,n are different , the DCP algorithm introduces some changes between the two cases because of the auxiliary term S 2. After a few simple steps, for Ampère law we find
where ε ‖,n=(1-d)ε 1+f[ε 1 n 2+ε ∞(1-n 2)] and g ‖,n=d(1-n 2), and for Faraday law
We want to replace E 2 with E in Eq. (3), but we cannot use directly Eq. (15) as in the case of Fig. 1(a). Since E 2 represents the field component along the integration direction of Faraday law, the projection technique  tells us that
By defining the auxiliary term S
In Fig. 1(c) we have a crossing only for the integration line of Faraday law so that in Eqs. (6)–(13) we set ε ‖,m=ε ∞ and c ‖,m=1. We thus obtain the standard FDTD algorithm where Eqs. (2) and (3) for E, D and S have coefficients
The last case that we discuss is shown in Fig. 1(f), corresponding to Fig. 1(b) with f=0 that implies ε ⊥,m=ε 1 and g ⊥,m=0. Because the coefficients in Eq. (20) contain the components of m, not used here, we have to recalculate the expression for E 2 using the components of n [see Eq. (17)]. Following the previous procedure leads to
The examples displayed in Fig. 1 are sufficiently general to show how to handle any further special case. Besides those already discussed, we just mention that a tilted flat interface would be treated in the same manner, except that n=m when two crossings occur. Actually, if the mesh is sufficiently fine, it is not a bad approximation to assume n=m also for a curved interface to simplify the implementation of Eqs. (13) and (20) [see Figs. 1(a) and 1(b)].
3. Numerical tests and discussion
We test the DCP against the standard staircase FDTD algorithm by computing the plasmon spectrum of a single metal NP  and that of two close-by metal NPs . These systems have the advantage of having an analytical solution in terms of vector spherical harmonics such that our results can be compared with the exact ones. Moreover, metal NPs have been shown to be a difficult task for FDTD [23, 24, 25] because their optical properties strongly depend on their shape and because staircasing introduces spurious peaks [25, 26].
Since we are working in two dimensions, we actually deal with a gold nanocylinder of radius r=25 nm placed in a background medium with index n b=1.7. The optical constants of gold  have been fitted to a Drude model with parameters ε ∞=9.9527, ω p=1.365×1016 rad/s and γ=1.2499×1014 rad/s. The fit becomes inaccurate for wavelengths shorter than about 600 nm , but that is not an issue as long as we are only interested in comparing FDTD with the exact solution, which is obtained using the same Drude model. The total scattering cross section (SCS) is computed using the total-field/scattered-field technique , as sketched in the inset of Fig. 2(c). A linearly polarized plane wave is incident with the magnetic field parallel to the nanocylinder axis. The total scattered power is obtained by summing the scattered power collected by a closed line in the scattered-field region. The FDTD mesh is terminated with PML absorbing boundary conditions  to ensure that the scattered field escapes without reflections at the mesh walls.
Figures 2(a) and 2(b) show the total SCS computed for a fine (Δ=0.5 nm) and a coarser (Δ=1.5 nm) mesh, respectively, using Mie theory , the staircase and the DCP FDTD. For a fine mesh both FDTD algorithms exhibit a very good agreement with the exact solution. For Δ=1.5 nm both methods are less accurate at the plasmon resonance, but staircasing shows an additional spurious peak at the wavelength λ≃625 nm. That does not only increase the discrepancy with the Mie results, but it may also lead to wrong conclusions, e.g., the existence of two plasmon resonances. The DCP algorithm cancels the spurious peak improving the accuracy with respect to the staircase FDTD.
The accuracy of an FDTD calculation becomes more critical if the metal losses decrease. To show that, we have computed the total SCS for halved γ in the Drude model [see Fig. 2(c)]. Notice that although the mesh pitch is only 0.4 nm, the plasmon spectrum has many more spurious peaks than in the previous case. The DCP is able to reduce the strength of such resonances yielding better agreement with the exact solution.
We study the performances of the DCP algorithm in more detail by looking at the relative error for the curves in Figs. 2(b) and 2(c). The results are plotted in Figs. 3(a) and 3(b), respectively. For Δ=1.5 nm and γ unchanged the strongest error is located at the spurious peak. The same occurs for Δ=0.4 nm and γ halved, but in this case there are very many oscillations and more spurious peaks. For both situations DCP exhibits an error which is less oscillatory and closer to zero than staircasing, except for a few situations. For instance DCP creates a spurious peak at λ=525 nm, where staircasing does not [see Fig. 3(b)].
To monitor the global performance of the two FDTD algorithms we have computed the average relative error over the wavelength range between 450 and 750 nm for several values of Δ. This error reads ∑λ|errλ|/N, where errλ is the wavelength-dependent relative error and N is the number of wavelength points in the sum . As shown in Fig. 3(c), the average error for DCP is clearly smaller than that for staircasing, especially for coarser meshes. The average error for the case of halved γ (dashed lines) is larger for both algorithms, but again DCP performs better than staircasing.
Besides losses, another factor that might significantly affect the accuracy of the FDTD algorithm is the presence of a strong and highly confined electric field such as in the case of two metal NPs separated by a small gap [22, 30, 31]. For this reason we compute the total SCS for a system composed of two gold nanocylinders with radius r=20 nm separated by a variable distance. The NPs are embedded in a dielectric background with index n b=1.7. A linearly polarized plane wave is incident with the magnetic field parallel to the nanocylinder axis and the wave vector oriented like the arrow sketched in the inset of Fig. 4(b).
Figure 4(a) displays the total SCS computed when the NPs have a center-to-center distance d=50 nm, meaning that the gap is 10 nm, and the mesh pitch is Δ=1.2 nm. The exact solution is given by generalized Mie theory implemented in the multiple multipole program (MMP) . While staircasing exhibits a strong spurious peak around λ=625 nm, DCP is always performing well except in between the two plasmon resonances, where staircasing reproduces the dip better. Figure 4(b) presents the same case for d=80 nm. Again, DCP is more accurate than staircasing and the agreement with the MMP calculation has improved for both methods. In fact for larger separations the near field is less confined and, consequently, its representation on the FDTD mesh as a step constant becomes more appropriate.
Like for a single NP we summarize the results by showing the average relative error computed for various discretizations [see Fig. 4(c)]. Notice that for d=50 nm the error in the DCP algorithm is comparable to staircasing, even though its oscillation as a function of Δ is smaller. On the other hand, for d=80 nm DCP becomes clearly better than staircasing, as found for a single NP. We interpret this fact by recalling that small separations generate such a strong and inhomogeneous field that its representation on the FDTD mesh as a step constant affects the accuracy as much as material interfaces. Furthermore, in Fig. 4(c) the average relative error exhibits a local minimum when the discretization is equal to 1.4 nm. We believe that this is a situation analogous to that of Fig. 6 in . However, the system we have here is too complicated to find out why this particular mesh pitch yields a smaller error.
We would like to emphasize that the DCP does not consider two partially filled cells that do not cut the integration line, i.e. the one that occupies a small space at the left bottom corner and the other one that almost fills the whole cell up to the right upper corner, see Fig. 1. By performing extensive numerical tests we have noticed that when the material interface has many of these cells, the DCP method gives the same results as staircasing. The chance for this to occur is small though, but it is higher for tiny objects such as NPs with a 10 nm diameter.
We have extended the CP-EP FDTD algorithm  to dispersive materials using the Z-transform formalism and adopting the Drude dispersion model as working case. We stress that the formulas in Eqs. (13)–(26) are valid for any linear dispersion model that requires terms up to z -2 S in Eq. (3). Furthermore, the DCP algorithm can be extended with little effort to treat nonlinear dispersive materials . Our method preserves the computational efficiency of the standard FDTD and it can be implemented by simply adding a pre-processing part to an existing FDTD code, without acting on the core subroutines.
We have compared the DCP algorithm with the staircase FDTD and the exact solutions for single gold NPs and NP pairs. We have shown that our method performs better than staircasing except in situations where the inaccuracy due to material interfaces is comparable with other sources of error, such as intense and strongly inhomogeneous near fields, or when many partially filled cells are neglected by the DCP algorithm. In these cases a possible solution would be subgridding . Besides an overall better accuracy with respect to staircasing, one of the major advantages of our approach is that it removes spurious peaks from plasmon spectra, which cause large discrepancies even for relatively fine meshes. The DCP algorithm is thus a more reliable FDTD method for studying the electromagnetic properties of complex metallo-dielectric nanostructures.
We thank Vahid Sandoghdar for continuous support and encouragement and Christian Hafner for helping with the MMP simulations. A. Mohammadi is thankful to the Persian Gulf University Research Council for partial support of this work. This work was also funded by the ETH Zurich initiative on Composite Doped Metamaterials (CDM).
References and links
1. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, 2005).
2. J.-Y. Lee and N.-H. Myung, “Locally tensor conformal FDTD method for modeling arbitrary dielectric surfaces,” Microwave Opt. Technol. Lett. 23, 245–249 (1999). [CrossRef]
3. A. Ditkowski, K. Dridi, and J. S. Hesthaven, “Convergent Cartesian grid methods for Maxwell’s equations in complex geometries,” J. Comput. Phys. 170, 39–80 (2001). [CrossRef]
4. W. Yu and R. Mittra, “A conformal finite difference time domain technique for modeling curved dielectric surfaces,” Microwave Opt. Technol. Lett. 11, 25–27 (2001).
5. K. H. Dridi, J. S. Hesthaven, and A. Ditkowski, “Staircase-free finite-difference time-domain formulation for general materials in complex geometries,” IEEE Trans. Antennas Propag. 49, 749–756 (2001). [CrossRef]
6. T. I. Kosmanis and T. D. Tsiboukis, “A systematic and topologically stable conformal finite-difference time-domain algorithm for modeling curved dielectric interfaces in three dimensions,” IEEE Trans. Microwave Theory Tech. 51, 839–847 (2003). [CrossRef]
7. J. Nadobny, D. Sullivan, W. Wlodarczyk, P. Deuflhard, and P. Wust, “A 3-D tensor FDTD-formulation for treatment of slopes interfaces in electrically inhomogeneous media,” IEEE Trans. Antennas Propag. 51, 1760–1770 (2003). [CrossRef]
8. A. Mohammadi, H. Nadgaran, and M. Agio, “Contour-path effective permittivities for the two-dimensional finite-difference time-domain method,” Opt. Express 13, 10367–10381 (2005) http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-25-10367. [CrossRef]
9. 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) http://www.opticsinfobase.org/abstract.cfm?URI=ol-31-20-2972. [CrossRef]
10. W. J. Padilla, D. N. Basov, and D. R. Smith, “Negative refractive index metamaterials,” Mat. Today 9, 28–35 (2006). [CrossRef]
11. V. M. Shalaev, “Optical negative-index metamaterials,” Nat. Photonics 1, 41–48 (2007). [CrossRef]
13. S. A. Maier and H. A. Atwater, “Plasmonics: localization and guiding of electromagnetic energy in metal/dielectric structures,” J. Appl. Phys. 98, 011101 (2005). [CrossRef]
15. R. Zia, J. A. Schuller, A. Chandran, and M. L. Brongersma, “Plasmonics: the next chip-scale technology,” Mat. Today 9, 20–27 (2006). [CrossRef]
16. A. Mohammadi and M. Agio, “Dispersive contour-path finite-difference time-domain algorithm for modeling surface plasmon polaritons at flat interfaces,” Opt. Express 14, 11330–11338 (2006) http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-23-11330. [CrossRef]
17. 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]
18. A. Deinega and I. Valuev, “Subpixel smoothing for conductive and dispersive media in the finite-difference time-domain method,” Opt. Lett. 32, 3429–3431 (2007) http://www.opticsinfobase.org/abstract.cfm?URI=ol-32-23-3429. [CrossRef]
19. D. M. Sullivan, “Z-transform theory and the FDTD method,” IEEE Trans. Antennas Propag. 44, 28–34 (1996). [CrossRef]
20. M. Born and E. Wold, Principles of Optics, 7th ed. (Cambridge U. Press, 1999).
21. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley Interscience, 1983).
22. H. Xu, J. Aizpurua, M. Käll, and P. Apell, “Electromagnetic contributions to single-molecule sensitivity in surface-enhanced Raman scattering,” Phys. Rev. E 62, 4318–4324 (2000). [CrossRef]
23. C. Oubre and P. Nordlander, “Optical properties of metallodielectric nanostructures calculated using the finite difference time domain method,” J. Phys. Chem. B 108, 17740–17747 (2004). [CrossRef]
24. C. Oubre and P. Nordlander, “Finite-difference time-domain studies of the optical properties of nanoshell dimers,” J. Phys. Chem. B 109, 10042–10051 (2005). [CrossRef]
25. F. Kaminski, V. Sandoghdar, and M. Agio, “Finite-difference time-domain modeling of decay rates in the near field of metal nanostructures,” J. Comput. Theor. Nanosci. 4, 635–643 (2007).
26. A. C. Cangellaris and D. B. 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–1524 (1991). [CrossRef]
27. CRC Handbook of Chemistry and Physics, 87th ed., D. R. Lide, ed. (CRC-Press, 2006) http://www.hbcpnetbase.com.
28. A. Vial, A.-S. Grimault, D. Macías, D. Biarchesi, and M. Lamy de la Chapelle, “Improved analytical fit of gold dispersion: application to the modeling of extinction spectra with a finite-difference time-domain method,” Phys. Rev. B 71, 085416 (2005). [CrossRef]
29. O. Ramadan and A. Y. Oztoprak, “Z-transform implementation of the perfectly matched layer for truncating FDTD domains,” IEEE Microwave Wirel. Compon. Lett. 13, 402–404 (2003). [CrossRef]
30. P. J. Schuck, D. P. Fromm, A. Sundaramurthy, G. S. Kino, and W. E. Moerner, “Improving the mismatch between light and nanoscale objects with gold bowtie nanoantennas,” Phys. Rev. Lett. 94, 017402 (2005). [CrossRef]
32. Ch. Hafner, Post-Modern Electromagnetics: Using Intelligent MaXwell Solvers, (John Wiley & Sons, 1999).