Effective Permittivity for FDTD Calculation of Plasmonic Materials

We present a new effective permittivity (EP) model to accurately calculate surface plasmons (SPs) using the finite-difference time-domain (FDTD) method. The computational representation of physical structures with curved interfaces causes inherent errors in FDTD calculations, especially when the numerical grid is coarse. Conventional EP models improve the errors, but they are not effective for SPs because the SP resonance condition determined by the original permittivity is changed by the interpolated EP values. We perform FDTD simulations using the proposed model for an infinitely-long silver cylinder and gold sphere, and the results are compared with Mie theory. Our model gives better accuracy than the conventional staircase and EP models for SPs.


Introduction
Surface plasmons (SPs) are waves that propagate along the surface of a metal and play a role in subwavelength imaging [1].SPs have been applied to advanced optical devices, for example, sensors, microscopies, lasers, invisibility cloaks, and hyperlenses [2][3][4][5][6][7][8][9].Numerical simulations are useful to analyze the complicated structures of these devices.In this paper, we discuss a finite-difference time-domain (FDTD) analysis of the SP.The FDTD method has gained popularity because of several reasons: it is easy to implement, it works in the time domain, and arbitrary shapes can be calculated [10][11][12][13].However, there are inherent errors in the FDTD method caused by numerical dispersion, the absorbing boundary condition, and the permittivity model.The former two can be improved by using the nonstandard FDTD method [14] and the perfectly matched layer [13].We address the error of the permittivity model.On a discrete grid, it is difficult to model permittivity interfaces that run between the grid points.The permittivity interface model is often an accuracy bottleneck, especially when the SPs are calculated on a coarse grid.By improving the model, we can accurately calculate optical effects on a coarse grid with a low computational cost.
The permittivity model error has been reduced by various approaches.The conformal FDTD method on non-orthogonal grids [15][16][17], counter path FDTD method [18,19], and subgridding method [20][21][22] can represent curved interfaces suitably, but they are relatively difficult to implement and increase the memory and computation time.A different approach using effective permittivities (EPs), which derives from interface interpolations based on Ampere's and Faraday's integration laws, can reduce the error of the permittivity model on coarse grids in a simple implementation and at low computational cost [23][24][25][26][27][28][29][30][31][32][33][34][35].The EP method, however, is not effective for plasmonic materials because the SP resonance condition is changed on interfaces due to interpolated values of EPs.Plasmonic waveguides have been analyzed by the conformal FDTD method using EPs [36], but the SP resonance wavelength is shifted several nanometers even on a fine grid (wavelength/grid spacing = 60).Furthermore, this method involves a complicated implementation given by a fourth-order differential equation.Ultimately, permittivities must be staircased to capture the correct SP resonance condition.
In this paper we propose a modified EP model to accurately calculate SP resonances on a coarse grid.A key concept is that the permittivities of plasmonic materials are used not only inside the materials but also on grid points adjacent to the materials (discussed in Section 2).In Section 3, numerical simulations for an infinitely-long silver cylinder and a gold sphere demonstrate that our model is in better agreement with Mie theory than the conventional staircase and EP models.

Effective Permittivity for Plasmonic Materials
Permittivities of plasmonic materials must be staircased, because the SP resonance condition given by interpolated EPs on interfaces is different from that of the original metal permittivity.But the staircase model is too crude to give high accuracy.Regardless of how close it lies to the interface, a grid point is either "in" or "out".The key to improve the representation is to take distance from the interface into account when assigning an permittivity value to a grid point.
In order to introduce our model, we start by analyzing EPs for plasmonic materials in the two-dimensional transverse magnetic (TM) mode, in which the electromagnetic fields reduce to three non-zero components E x , E y , and H z .Ampere's and Faraday's integration laws are given by ∂ ∂t where E is the electric field, D is the electric induction, H is the magnetic field, B is the magnetic induction, dS is an element of surface area, and dl is an infinitesimal length element.Figure 1 shows examples of partially filled FDTD cells in the TM mode: (a) interface parallel to E x ; (b) interface perpendicular to E x , where ∆x, ∆y are grid spacings, d, f are filling factors on the integration lines, and ε 1 , ε 2 are permittivities outside and inside the plasmonic material, respectively.In Figure 1(a), Ampere's law gives an integration with respect to the E x components: Assuming ∆y to be sufficiently small, D x and E x are essentially constant in the integration region.Since E x parallel to the interface is continuous on each plane, substituting D x = ε ∥ E x (ε ∥ = EP for Ampere's law) into the left-hand side of Equation ( 3), we obtain In Figure 1(b), Faraday's law gives the other integration: Since the D x perpendicular to the interface is continuous, similarly we can obtain the EP for Faraday's law: Detailed derivations of ε ∥ and ε ⊥ are given in [23][24][25][26][27][28].
The ε ∥ and ε ⊥ models give high accuracy for non-plasmonic materials on a coarse grid, however, they are not effective for plasmonic ones.Figure 2 shows (a) propagating and (b) localized SPs.In Figure 2, strongly resonant SP fields enclosed with dotted lines are in the polarization perpendicular to the interface.The ε ⊥ model fits this situation, but the approximation of Equation ( 6) is invalid because the SP decays exponentially with increasing distance from the metal surface.On the other hand, the damping has only a small effect for the low power parallel field, and the ε ∥ model can be used.However, the permittivity values must be staircased in order to preserve the resonance condition as discussed before.We empirically found that the SP resonances can be accurately calculated by staircasing the ε ∥ model in the following manner: for both the parallel and perpendicular cases.Equation ( 8) represents what we call the staircased EP (S-EP) model.In the S-EP model, even when the E x grid point is outside the plasmonic medium (ε 2 ) as shown in Figure 1(a), it is assigned to the value ε 2 if d > 0. When the E x is perpendicular to the interface as shown in Figure 1(b), the conventional staircase model is used because d = 0 or 1. Equation ( 8) is applied to the E x grid points, but the E y assignments can be obtained with the same procedure.The S-EP model can be simply implemented in Yee Algorithm [10][11][12][13]: where ∆t is the time interval, h is the grid spacing, µ 0 is the permeability, and we simply write D x (t = n∆t) → D n x (n = integer).The spatial difference operators d x , d y are defined by The S-EP model thus prescribes how to assign permittivity values to grid points near medium interfaces.In Section 3, we shown that the S-EP model is more accurate than the conventional staircase and EP models.This method can be readily extended to three dimensions, keeping in mind that the line integrals become surface integrals in Ampere's law.

Numerical Validations
We validate the S-EP model in FDTD simulations for an infinitely-long cylinder of silver (Ag) and sphere of gold (Au).Localized SPs are generated for wavelengths: >340 nm for the Ag cylinder; >510 nm for the Au sphere.FDTD results using the S-EP model are compared with those of the conventional staircase and EP models, and with Mie theory [37].Since the ε ⊥ model has a large error for SPs as discussed in Section 2, we employ the ε ∥ model as a comparison.To analyze dispersive materials, we use the auxiliary differential equation method with the Drude Model (taking an example) [13]: where ω p is the plasma frequency, and ν c is the collision frequency.To capture the measured permittivities of Johnson and Christy [38], we calculate ω p , ν c at each wavelength as follows: where ε = ε r + iε i .

Infinitely-Long Ag Cylinder
Figure 3 shows the simulation setup for an infinitely-long Ag cylinder: the computational domain is terminated by a perfectly matched layer (PML) [13]; a plane wave source of wavelength λ 0 = 430.501nm (696 THz) in the SP waveband of Ag (>340 nm) is TM polarized; the cylinder radius is a = 1.25λ 0 = 538.126nm; the grid spacing is h = 20 nm (λ 0 /h = 21.5); and the time step is given by the Courant limit.FDTD calculations are compared with Mie theory on a contour C of radius 1.2a about the cylinder.Simulation parameters for the infinitely-long Ag cylinder are listed in Table 1.We calculated the scattered intensity in the near field using the FDTD method until the steady state (about 50 wave periods).Figure 4    Figure 5 shows root-mean-square (RMS) errors of FDTD results using staircase, EP, and S-EP models relative to Mie theory on C. Figure 5(a) shows the RMS errors as a function of incident wavelength (h = 20 nm, a = 538.126nm).In the SP waveband (>340 nm) the accuracy of FDTD calculations is relatively low, but the S-EP model provides much higher accuracy than the staircase and EP models.As the grid spacing decreases the model errors converge, but the S-EP model is more accurate than the conventional staircase and EP models on a coarse grid.This greater accuracy on a coarse grid is the advantage of our method, because computational cost rises rapidly with decreasing grid size.Figure 5(c) shows RMS errors as a function of cylinder radius (h = 20 nm, λ 0 = 430.501nm).The accuracy of the S-EP model is independent of the structure size.This robustness against the change of material parameters indicates that the S-EP model is not just an optimization for a specific example.Finally, we validate the S-EP model for subwavelength structures.Figure 5(d) shows RMS errors as a function of subwavelength cylinder radius using a finer grid spacing (h = 10 nm, λ 0 = 430.501nm).The accuracy of the EP model depends strongly on the subwavelength structure size, but the S-EP model gives robust and accurate results.

Au Sphere
In three dimensions the EP and S-EP methods are given by exchanging line integrals for surface integrals in Ampere's law. Figure 6 shows the simulation setup for a Au sphere: the computational domain is terminated by the PML; a plane wave source of wavelength λ 0 = 616.837nm (367 THz) in the SP waveband of Au (>510 nm) is polarized along the y axis and propagates in the +z direction; the cylinder radius is a = 1.5λ 0 = 925.255nm; the grid spacing is h = 20 nm (λ 0 /h = 30.8);and the time step is given by the Courant limit.FDTD calculations are compared with Mie theory on a surface S of radius 1.2a and on a contour C on S in the y-z plane.Simulation parameters for the Au sphere are listed in Table 2.

Conclusions
We have proposed a S-EP model for FDTD computation of plasmonic materials.The S-EP model is a modification of the EP model derived from Ampere's law.FDTD simulations using the S-EP model were implemented for the infinitely-long Ag cylinder and Au sphere.Numerical results were compared with Mie theory in the SP waveband, and the S-EP model gave higher accuracy than the conventional staircase and EP models even on a coarse grid.The S-EP model are robust against changes of material parameters, which indicates that the methodology is valid over a broad spectrum of problems.
The S-EP model can be applied to other dispersion models and simulation methods, because it is implemented only by changing how permittivity values are assigned to grid points.We have obtained the same performance using the Lorentz model [13], the Drude model with critical points [39,40], the recursive convolution (RC) FDTD method [41], and the piecewise linear RC-FDTD method [42].

Figure 1 .
Figure 1.Partially filled FDTD cells in the TM mode: (a) interface parallel to E x ; (b) interface perpendicular to E x .Dotted lines represent the integration lines in Ampere's and Faraday's laws.d, f represent filling factors on the integration lines.ε 1 , ε 2 are permittivities outside and inside the plasmonic material, respectively.

Figure 2 .
Figure 2. Surface plasmon (SP) resonances: (a) propagating SP; (b) localized SP.Electric charge distributions and lines of force are depicted.Regions enclosed with dotted lines indicate strongly resonant SP fields.

Figure 3 .
Figure 3. Simulation setup for an infinitely-long Ag cylinder.Incident TM polarization (a = radius, h = grid spacing).Ag cylinder (gray region) is immersed in air (white region).FDTD results are compared with Mie theory on dotted contour C.
shows the calculation results for the E y scattered intensity: (a) Mie theory; (b) staircase model; (c) EP model; (d) S-EP model; and (e) angular distributions on C. As shown in Figure4, the calculated result using the S-EP model is in much better agreement with Mie theory than those of staircase and EP models.

Figure 4 .
Figure 4. Calculated E y scattered intensity distributions for an infinitely-long Ag cylinder using: (a) Mie theory; (b) Staircase model; (c) EP model; (d) S-EP model (black region represents the cylinder); (e) Angular intensity distributions on contour C of radius 1.2a.

Figure 5 (
Figure5shows root-mean-square (RMS) errors of FDTD results using staircase, EP, and S-EP models relative to Mie theory on C. Figure5(a) shows the RMS errors as a function of incident wavelength (h = 20 nm, a = 538.126nm).In the SP waveband (>340 nm) the accuracy of FDTD calculations is relatively low, but the S-EP model provides much higher accuracy than the staircase and EP models.Figure 5(b) shows the RMS errors as a function of grid spacing (λ 0 = 430.501nm, a = 538.126nm).As the grid spacing decreases the model errors converge, but the S-EP model is more accurate than the conventional staircase and EP models on a coarse grid.This greater accuracy on a coarse grid is the advantage of our method, because computational cost rises rapidly with decreasing grid size.Figure5(c)shows RMS errors as a function of cylinder radius (h = 20 nm, λ 0 = 430.501nm).The accuracy of the S-EP model is independent of the structure size.This robustness against the change of material parameters indicates that the S-EP model is not just an optimization for a specific example.Finally, we validate the S-EP model for subwavelength structures.Figure5(d)shows RMS errors as a function of subwavelength cylinder radius using a finer grid spacing (h = 10 nm, λ 0 = 430.501nm).The accuracy of the EP model depends strongly on the subwavelength structure size, but the S-EP model gives robust and accurate results.

Figure 6 .
Figure 6.Simulation setup for a Au sphere (a = radius, h = grid spacing).Incident E y polarized wave propagates along the z axis.Au sphere (gray region) is immersed in air (white region).FDTD results are compared with Mie theory on surface S of radius 1.2a and contour C on S in the y-z plane (dotted line).

Figure 7
Figure 7 shows calculated E y scattered intensities on S at the steady state: (a) Mie theory; (b) staircase model; (c) EP model; (d) S-EP model; and (e) angular distributions on C. As shown in Figure 7, the S-EP model gives close agreement with Mie theory better than the staircase and EP models.Consequently, the S-EP model is also effective for the three-dimensional SP resonance.

Figure 7 .
Figure 7. Calculated E y scattered intensity distributions for a Au sphere on surface S of radius 1.2a using: (a) Mie theory; (b) Staircase model; (c) EP model; (d) S-EP model; (e) Angular intensity distributions on contour C in the y-z plane.

Table 1 .
Simulation parameters for an infinitely-long Ag cylinder.

Table 2 .
Simulation parameters for a Au sphere.