## Dispersive contour-path finite-difference time-domain algorithm for modelling surface plasmon polaritons at flat interfaces

Optics Express, Vol. 14, Issue 23, pp. 11330-11338 (2006)

http://dx.doi.org/10.1364/OE.14.011330

Acrobat PDF (117 KB)

### Abstract

We investigate the accuracy of the two-dimensional Finite-Difference Time-Domain (FDTD) method in modelling Surface Plasmon Polaritons (SPPs) in the case of a single metal-dielectric interface and of a thin metal film showing that FDTD has difficulties in the low-group-velocity region of the SPP. We combine a contour-path approach with Z transform to handle both the electromagnetic boundary conditions at the interface and the negative dispersive dielectric function of the metal. The relative error is thus significantly reduced.

© 2006 Optical Society of America

## 1. Introduction

2. W.L. Barnes, A. Dereux, and T.W. Ebbesen, “Surface plasmon subwavelength optics,” Nature **424**, 824–830 (2003). [CrossRef] [PubMed]

4. D. Sullivan, *Electromagnetic Simulation Using the FDTD Method* (IEEE Press, 2000). [CrossRef]

5. T.G. Jurgens, A. Taflove, K. Umaschankar, and T.G. Moore, “Finite-Difference Time-Domain Modeling of Curved Surfaces,” IEEE Trans. Antennas Propag. **40**, 357–365 (1992). [CrossRef]

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

5. T.G. Jurgens, A. Taflove, K. Umaschankar, and T.G. Moore, “Finite-Difference Time-Domain Modeling of Curved Surfaces,” IEEE Trans. Antennas Propag. **40**, 357–365 (1992). [CrossRef]

8. K.-P. Hwang and A.C. Cangellaris, “Effective Permittivities for Second-Order Accurate FDTD Equations at Dielectric Interfaces,” IEEE Microwave Wirel. Compon. Lett. **11**, 158–160 (2001). [CrossRef]

9. A. Mohammadi and M. Agio, “Contour-path effective permittivities for the twodimensional finite-difference time-domain method,” Opt. Express **13**, 10367–10381 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-25-10367, and references therein. [CrossRef] [PubMed]

10. M.K. Kärkkäinen, Subcell FDTD Modeling of Electrically Thin Dispersive Layers,” IEEE Trans. Microwave Theory Tech. **51**, 1774–1780 (2003). [CrossRef]

11. D. Popovic and M. Okoniewski, “Effective Permittivity at the Interface of Dispersive Dielectrics in FDTD,” IEEE Microwave Wirel. Compon. Lett. **13**, 265–267 (2003). [CrossRef]

*staircase*approximation. Moreover, the field components must be displaced in space and time to construct the finitedifference equations [3]. For these reasons, the longitudinal and transverse electric components of a SPP will sense different interfaces. Therefore, regardless of the position, FDTD is unable to fit the correct conditions of a SPP even for a flat interface.

5. T.G. Jurgens, A. Taflove, K. Umaschankar, and T.G. Moore, “Finite-Difference Time-Domain Modeling of Curved Surfaces,” IEEE Trans. Antennas Propag. **40**, 357–365 (1992). [CrossRef]

9. A. Mohammadi and M. Agio, “Contour-path effective permittivities for the twodimensional finite-difference time-domain method,” Opt. Express **13**, 10367–10381 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-25-10367, and references therein. [CrossRef] [PubMed]

## 2. Dispersive contour-path method

*ε*

_{d}and

*ε*

_{m}(

*ω*)=

*ε*

_{∞}-

*ω*(

*ω*+

*ιη*)), respectively;

*ω*

_{p}is the plasma frequency and

*η*is the damping frequency [12]. The constitutive relations

**B**=

**H**and

**D**=

*ε*

_{m}(

*ω*)

**E**, where

*μ*

_{0}and

*ε*

_{0}are set to 1, must be expressed in time domain to be compatible with the FDTD paradigm. Therefore, the relationship between

**D**and

**E**becomes a convolution. We use the Z transform as a suitable transformation to replace discrete convolutions with multiplications [4

4. D. Sullivan, *Electromagnetic Simulation Using the FDTD Method* (IEEE Press, 2000). [CrossRef]

**D**and

**E**becomes for each field component

*t*is the time step used in FDTD and

*ε*

_{m}(

*z*) corresponds to

*ε*

_{m}(

*ω*) in z-space

4. D. Sullivan, *Electromagnetic Simulation Using the FDTD Method* (IEEE Press, 2000). [CrossRef]

*ε*

_{m}(

*ω*), i.e. by the dispersion model. Since in the Z space,

*z*

^{-1}works as a delay operator [4

4. D. Sullivan, *Electromagnetic Simulation Using the FDTD Method* (IEEE Press, 2000). [CrossRef]

^{n+1}represents the field at time

*t*=(

*n*+1)Δ

*t*.

9. A. Mohammadi and M. Agio, “Contour-path effective permittivities for the twodimensional finite-difference time-domain method,” Opt. Express **13**, 10367–10381 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-25-10367, and references therein. [CrossRef] [PubMed]

*E*

_{x},

*E*

_{y},

*H*

_{z}). Consider first the

*E*

_{y}component, which is parallel to the interface. The FDTD integration path associated to Ampére law is partially inside the metal and partially in the dielectric. For this reason, we introduce an average dielectric displacement

*d*and Δ are defined in Fig. 1,

*D*

_{m}and

*D*

_{d}refer to the electric displacement in metal and dielectric, respectively, and the

*y*subscript is omitted. Using Eq. (3) and exploiting the boundary condition

*E*

_{m}=

*E*

_{d}≡

*E*, we obtain

8. K.-P. Hwang and A.C. Cangellaris, “Effective Permittivities for Second-Order Accurate FDTD Equations at Dielectric Interfaces,” IEEE Microwave Wirel. Compon. Lett. **11**, 158–160 (2001). [CrossRef]

**13**, 10367–10381 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-25-10367, and references therein. [CrossRef] [PubMed]

*D*and

*E*

_{x}field component, which is perpendicular to the interface. The integration path of Faraday law is again partially inside the metal and partially in the dielectric, as shown in Fig. 1. Following the same idea, we introduce an average electric field

*f*is defined in Fig. 1 and the

*x*subscript is omitted. Using Eq. (3) and the boundary condition

*D*

_{m}=

*D*

_{d}≡

*D*, we obtain

8. K.-P. Hwang and A.C. Cangellaris, “Effective Permittivities for Second-Order Accurate FDTD Equations at Dielectric Interfaces,” IEEE Microwave Wirel. Compon. Lett. **11**, 158–160 (2001). [CrossRef]

**13**, 10367–10381 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-25-10367, and references therein. [CrossRef] [PubMed]

*S*

_{m}(

*z*), the electric field inside the metal is needed, which can be obtained from

*E*. After a few simple algebraic operations and by saving

*E*and

## 3. Numerical tests and discussion

*ω*(

*k*) of SPPs propagating along a single metal-dielectric interface and along a thin metal film. These systems are chosen for two reasons; first because they have an analytical solution, which we can take as reference [1]; second because they represent the simplest cases for understanding the basic issues of FDTD in handling SPPs.

### 3.1. Single metal-dielectric interface

13. C.T. Chan, Q.L. Yu, and K.M. Ho, “Order-*N* spectral method for electromagnetic waves,” Phys. Rev. B **51**, 16635–16642 (1995). [CrossRef]

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

13. C.T. Chan, Q.L. Yu, and K.M. Ho, “Order-*N* spectral method for electromagnetic waves,” Phys. Rev. B **51**, 16635–16642 (1995). [CrossRef]

*ε*

_{d}=2.25) and copper (

*ε*

_{∞}=1.0,

*ω*

_{p}=5.0×10

^{15}rad/s and

*η*=0.01

*ω*

_{p}) [15]. The dispersion relations computed using staircasing, the CP method and the analytical solutions are plotted in Fig. 2. While in the region where the SPP has a high group velocity both methods agree very well with the analytical curve, as the group velocity decreases staircasing exhibits a much larger error. Indeed, this is the case where the SPP is tightly localised at the interface and thus a proper treatment of the electromagnetic boundary conditions is more important. In fact, the relative error |

*ω*(

*k*) FDTD-

*ω*(

*k*)|/

*ω*(

*k*) increases linearly with

*k*for frequencies in the flat region of the dispersion relation (

*k*≥

*k*

_{s}), for both staircase and CP.

*N*is the number of wavevectors used in the sum [16] and

*ω*(

*k*) is the solution to Eq. (13). While the CP error remains small for a wide range of cell sizes, the staircasing one grows rapidly and becomes larger than 10% already for discretizations as small as Δ=10nm. By fitting the curves of Fig. 3 with the function

*y*=

*a*Δ

^{α}we can estimate the order of convergence [17]. While staircasing is only first-order accurate (

*α*=0.95±0.02), the CP-FDTD exhibits

*α*=1.47±0.1. We have found basically the same values for

*α*by fitting the

*L*

^{2}error

*α*=2) as it occurs for interfaces between dielectric materials [8

**11**, 158–160 (2001). [CrossRef]

*η*=0, as trial solution we take the following SPP mode

*i*,

*j*point to an FDTD cell. Notice that

*k*is exact since it is set by the Bloch boundary conditions, while the lateral confinement

*k̃*and the frequency

*are affected by the discretization Δ. Moreover,*ω ˜

*k̃*and

*E*

_{x,0}are different in the dielectric and in the metal. After substituting Eqs. (14) into the Yee algorithm, eliminating the common exponential factors and applying the boundary conditions at the metal-dielectric interface, we obtain the expression for the SPP dispersion relation on the FDTD mesh,

*S*

_{c}is the Courant stability factor [3], which in our numerical experiments is 0.95. By computing the relative error for different discretizations Δ and fitting it with

*y*=

*a*Δ

^{α}, as we have done for the data in Fig. 3, we obtain

*α*=2, like for a homogeneous medium. In deriving Eq. (15) we assumed that the boundary conditions are fulfilled and that the dielectric function in the metal is exactly

*ε*

_{m}(

*) at*ω ˜

*. Since material dispersion in FDTD enters through a convolution between*ω ˜

*ε*

_{m}(

*n*Δ

*t*) and the electric field, which is also affected by the discretization, we expect also an error for

*ε*(

*) at*ω ˜

*. Following the guidelines of Ref. [18*ω ˜

18. C. Hulse and A. Knoesen, “Dispersive models for the finite-difference time-domain method: design, analysis, and implementation,” J. Opt. Soc. Am. A **11**, 1802–1811 (1994). [CrossRef]

*t*

^{2}. Therefore, Eq. (15) should be corrected as

*y*=

*a*Δ

_{α}to estimate the order of convergence.We find that

*α*remains equal to 2, which makes sense given that the error on

*ε*(

*ω*) is second-order in Δ

*t*. The fact that

*α*does not reach 2 for the CP method should therefore be attributed to other sources of error. We very much believe that the major problem resides in the Bloch boundary conditions, which are not fully correct for a system with losses. Indeed, by testing our method on a system with significantly lower losses we obtained

*α*≃1.8. We could not test a system without losses because FDTD must fulfil the causality condition for the dielectric function. Despite this problem we are confident that the relative improvement of the CP method against staircasing remains valid, because the boundary conditions affect both algorithms equally.

*d*and

*f*shown in Fig. 1. Actually,

*d*and

*f*are not independent, because

*d*=

*f*+0.5Δ if

*f*<0.5Δ or

*d*=

*f*-0.5Δ if

*f*≥0.5Δ. Therefore, it is enough to study the case when

*f*varies from 0 to Δ. Notice that staircasing implies only two situations

*f*<0.5Δ and

*f*≥0.5Δ that exhibit different errors as it can be inferred from Figs. 2 and 3. On the contrary, the CP method is able to account for the exact value of

*f*. This property is very important for optimisation algorithms where a structure parameter is finely tuned. The result, plotted in Fig. 4, shows that the error depends on

*f*in a non-trivial manner. The error is larger when

*f*/Δ is about 0.75 and, like for staircasing, it is positive for

*f*/Δ≃0.25 and negative for

*f*/Δ≃0.75. The minimum error is found for

*f*/Δ≃0.5, corresponding to the situation where the cell associated to the parallel field component is completely filled with metal, and for

*f*/Δ≃1.0 or 0.0, corresponding to the situation where the cell associated to the perpendicular field component is completely filled with metal. In other words, the best performances are obtained when the dispersive CP operates only on one field component. In this perspective, one could think of positioning the mesh so that only one CP integral is used, i.e. only one field component has a partially filled cell [11

11. D. Popovic and M. Okoniewski, “Effective Permittivity at the Interface of Dispersive Dielectrics in FDTD,” IEEE Microwave Wirel. Compon. Lett. **13**, 265–267 (2003). [CrossRef]

*f*, we believe that avoiding such constraint would be more advantageous for modelling general structures.

### 3.2. Thin metal film with symmetric interfaces

*d*=50nm) in glass. With respect to the previous case we have two new features: the interaction of SPPs bound at the two interfaces and the existence of a dimension for the structure i.e. the film thickness. The staircasing error therefore arises both from the individual interfaces and from the imperfect matching of the film thickness, which is responsible for the SPP coupling.

*ω*

_{+}) and symmetric (

*ω*

_{-}) SPPs. Their dispersion relations are given in implicit form by the following equations [1],

*f*=0.75Δ, the lower one will have

*f*=0.25Δ, and vice versa. Because we see from Fig. 2 that one curve is above and the other one is below the correct result, it could happen that their coupling into symmetric and antisymmetric modes luckily reduces the error, explaining the dip in Fig. 6. The same occurs for the CP method since it also exhibits the change of sign in the error as shown in Fig. 4. Notice that this effect occurs also for Δ=2nm, which shows up as a steeper drop in the error. It is not so evident for the CP method though, because for such discretization we have already reached the error on the Fourier transform. These findings suggest that under certain conditions one could reduce the error more rapidly than at the rate set by

*α*. On the other hand, such route seems not realistic for complicated structures, where the various cancelling effects cannot be guided by intuition. What is important here is that, again, the CP approach is always nearly an order of magnitude better than staircasing.

## 4. Conclusions

## Acknowledgements

## References and links

1. | H. Raether, |

2. | W.L. Barnes, A. Dereux, and T.W. Ebbesen, “Surface plasmon subwavelength optics,” Nature |

3. | A. Taflove and S.C. Hagness, |

4. | D. Sullivan, |

5. | T.G. Jurgens, A. Taflove, K. Umaschankar, and T.G. Moore, “Finite-Difference Time-Domain Modeling of Curved Surfaces,” IEEE Trans. Antennas Propag. |

6. | C. Oubre and P. Nordlander, “Optical Properties of Metallodielectric Nanostructures Calculated Using the Finite Difference Time Domain Method,” J. Phys. Chem. B |

7. | H. Shin and S. Fan, “All-Angle Negative Refraction for Surface Plasmon Waves Using a Metal-Dielectric-Metal Structure,” Phys. Rev. Lett. |

8. | K.-P. Hwang and A.C. Cangellaris, “Effective Permittivities for Second-Order Accurate FDTD Equations at Dielectric Interfaces,” IEEE Microwave Wirel. Compon. Lett. |

9. | A. Mohammadi and M. Agio, “Contour-path effective permittivities for the twodimensional finite-difference time-domain method,” Opt. Express |

10. | M.K. Kärkkäinen, Subcell FDTD Modeling of Electrically Thin Dispersive Layers,” IEEE Trans. Microwave Theory Tech. |

11. | D. Popovic and M. Okoniewski, “Effective Permittivity at the Interface of Dispersive Dielectrics in FDTD,” IEEE Microwave Wirel. Compon. Lett. |

12. | M. Born and E. Wolf, |

13. | C.T. Chan, Q.L. Yu, and K.M. Ho, “Order- |

14. | O. Ramadan and A.Y. Oztoprak, “Z-transform implementation of the perfectly matched layer for truncating FDTD domains,” IEEE Microwave Wirel. Compon. Lett. |

15. | C. Kittel, |

16. |
In the error calculation we choose wavevectors |

17. |
The fit is performed using the logarithm of the error, so that the exponent |

18. | C. Hulse and A. Knoesen, “Dispersive models for the finite-difference time-domain method: design, analysis, and implementation,” J. Opt. Soc. Am. A |

**OCIS Codes**

(000.4430) General : Numerical approximation and analysis

(130.0130) Integrated optics : Integrated optics

(240.6680) Optics at surfaces : Surface plasmons

**ToC Category:**

Optics at Surfaces

**History**

Original Manuscript: September 5, 2006

Revised Manuscript: October 25, 2006

Manuscript Accepted: October 25, 2006

Published: November 13, 2006

**Citation**

Ahmad Mohammadi and Mario 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/oe/abstract.cfm?URI=oe-14-23-11330

Sort: Year | Journal | Reset

### References

- H. Raether, Surface Plasmons on Smooth and Rough Surfaces and on Gratings (Springer-Verlag, 1988)
- W. L. Barnes, A. Dereux, and T. W. Ebbesen, "Surface plasmon subwavelength optics," Nature 424,824-830 (2003). [CrossRef] [PubMed]
- A. Taflove, and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, Norwood, MA, 2005).
- D. Sullivan, Electromagnetic Simulation Using the FDTD Method (IEEE Press, 2000). [CrossRef]
- T. G. Jurgens, A. Taflove, K. Umaschankar, and T. G. Moore, "Finite-difference time-domain modeling of curved surfaces," IEEE Trans. Antennas Propag. 40,357-365 (1992). [CrossRef]
- 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]
- H. Shin, and S. Fan, "All-angle negative refraction for Surface Plasmon Waves using a metal-dielectric-metal structure," Phys. Rev. Lett. 96, 073907(4) (2006).
- K.-P. Hwang, and A. C. Cangellaris, "Effective permittivities for second-order accurate FDTD equations at Dielectric interfaces," IEEE Microw. Wirel. Compon. Lett. 11,158-160 (2001). [CrossRef]
- A. Mohammadi, and M. Agio, "Contour-path effective permittivities for the twodimensional finite-difference time-domain method," Opt. Express 13,10367-10381 (2005). [CrossRef] [PubMed]
- M. K. Kärkkäinen, Subcell FDTD modeling of electrically thin dispersive layers," IEEE Trans. Microwave Theory Tech. 51,1774-1780 (2003). [CrossRef]
- D. Popovic, and M. Okoniewski, "Effective permittivity at the interface of dispersive dielectrics in FDTD," IEEE Microw. Wirel. Compon. Lett. 13,265-267 (2003). [CrossRef]
- M. Born and E. Wolf, Principles of Optics (Cambridge U. Press, Cambridge, UK, 1999), 7th ed.
- C. T. Chan, Q. L. Yu, and K. M. Ho, "Order-N spectral method for electromagnetic waves," Phys. Rev. B 51,16635-16642 (1995). [CrossRef]
- O. Ramadan and A. Y. Oztoprak, "Z-transform implementation of the perfectly matched layer for truncating FDTD domains," IEEE Microw. Wirel. Compon. Lett. 13,402-404 (2003). [CrossRef]
- C. Kittel, Introduction to Solid State Physics (Wiley, New York, 1996), 7th ed.
- In the error calculation we choose wavevectors k ≥ 1.5ks (Fig. 3 and Fig. 4) and k ≥ 2.5ks (Fig. 6) to focus on the low-group-velocity region.
- The fit is performed using the logarithm of the error, so that the exponent 〈 is obtained from log y = 〈 log x+loga. The errors for ⊗ = 2 −4nm are excluded from the fit since the Fourier error is of the same order of magnitude and thus it may have an effect on the actual FDTD accuracy.
- C. Hulse and A. Knoesen, "Dispersive models for the finite-difference time-domain method: design, analysis, and implementation," J. Opt. Soc. Am. A 11, 1802-1811 (1994). [CrossRef]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.