## Abstract

We introduce a simple scattered field (SF) technique that enables finite difference time domain (FDTD) modeling of light scattering from dispersive objects residing in stratified dispersive media. The introduced SF technique is verified against the total field scattered field (TFSF) technique. As an application example, we study surface plasmon polariton enhanced light transmission through a 100nm wide slit in a silver film.

©2010 Optical Society of America

## 1. Introduction

The finite difference time domain (FDTD) method [1] is widely used as a numerical tool to solve electromagnetic scattering problems, where a spatially finite scatterer is illuminated by a plane wave. Plane waves are typically excited in FDTD either by the total field scattered field (TFSF) [2] or by the scattered field (SF) technique [3]. In the absence of a scatterer, both of these techniques produce a field that is called the incident field. When a scatterer interacts with the incident field, a scattered field is formed. The sum of the incident and the scattered field is referred to as the total field. The TFSF technique divides the computation domain into the total field and the scattered field region separated by an interface that generates the incident field and guides it through the total field region. In the SF technique, Maxwell’s equations are reformulated in such a manner that only the scattered fields are calculated by the FDTD method. Since the incident field is initially known, the total field is obtained as a sum of the incident and the scattered field.

Excitation of the incident field in FDTD is straightforward when scatterers reside in a homogenous background medium. The situation, however, is more complicated when the scatterers are located in a stratified background structure. This is mainly due to the fact that modeling propagation of the incident field in a uniform medium is much simpler than in stratified media. An exception is a case where a normally incident plane wave illuminates stratified media since an auxiliary one-dimensional FDTD simulation can be used as a real-time look-up table for the TFSF technique [2]. To handle obliquely incident plane waves in the TFSF technique, Winton et al. [4] have presented an approach based on multiple one-dimensional FDTD grids. The approach was further developed by Jian et. al. [5]. Schneider and Abdijalilov [6,7] have calculated incident fields for TFSF boundary analytically, including the numerical dispersion effects to obtain a non-leaking boundary. Demarest et al. [8] have introduced the SF technique for oblique plane waves in non-dispersive stratified media.

To our knowledge, this is the first time when the SF technique is presented for dispersive scatterers in dispersive stratified media. A SF formulation for dispersive materials was presented by Kong et al. [9], but in their formulation it is assumed that the incident field propagates in a uniform medium. In addition, their SF formulations depend on the used dispersive medium (Lorentz/Debye/Drude). In this paper, we present a simple approach that is nearly independent of the used dispersive medium.

This paper is organized as follows. In the next section, we introduce the scattered field technique in the frequency domain and then transform the obtained equations to the time domain. Section 3 considers in detail how the source term of the SF technique is calculated in practice. In section 4, we introduce the SF technique for scatterers in perfect electric conductor layers. The SF technique is then verified against the TFSF technique in Section 5. Surface plasmon polariton enhance light transmission through a slit is studied in Section 6 via the introduced SF technique. Finally, the summary and conclusions are presented in Section 7.

## 2. Scattered field technique

This section introduces the SF technique for dispersive scatterers in a dispersive stratified background structure. The SF technique is based on the fact that due to linearity of Maxwell’s equations, the total electric and the magnetic field can be expressed as a sum of the incident field and the scattered field. That is,

We start from the frequency domain Maxwell’s equations for a linear, non-magnetic (*μ* = *μ*_{0}) medium:

*μ*

_{0}is the permeability of free space, and

*ε*(

*ω*) =

*ε*

_{re}(

*ω*) +

*jε*

_{im}(

*ω*) is the complex dielectric permittivity. We assume that the incident fields obey the following equations:

*ε*

_{inc}(

*ω*) is the local complex dielectric permittivity of the structure for which the incident field is known. Substituting Eqs. (1) and (2) into Eqs. (3) and (4), using Eqs. (5) and (6), and then rearranging terms, we have

Maxwell’s equations are now reformulated to such a form that only the scattered field components are required to be solved with the assumption that the incident field is initially known. Note that the source term **S**(*ω*), which is the only quantity that depends on the incident field, is non-zero only in regions where *ε*(*ω*) differs from *ε*_{inc}(*ω*).

The frequency domain Eqs. (7)-(9) must be transformed to the time domain before they can be used in the FDTD method. Using the inverse Fourier transform, we obtain

*F*(

^{−1}{jωε*ω*)

**E**

_{sca}(

*ω*)

*}*is typically incorporated into a FDTD model via the auxiliary differential equation (ADE) [10] or the recursive convolution method [11]. The SF technique presented herein is compatible with the both approaches. A total field FDTD implementation can be easily converted to calculate only scattered fields. The only required change is that the ∇ ×

**H**(

*t*) term is replaced in the total electric field update equation by ∇ ×

**H**

_{sca}(

*t*) −

**(**

*S**t*).

## 3. Calculation of the source term

To calculate the source term (12), first the interaction of the illuminating plane wave with the stratified background structure has to be solved in the frequency domain. Throughout this paper, we assume that the illuminating plane wave pulse exhibits the following time dependency in the coordinate origin:

*f*

_{0}is the central frequency of the pulse,

*t*

_{0}is the time domain displacement of the pulse peak from the origin, and

*τ*defines the 1/

*e*half width of the pulse in the time domain. First,

*f*(

*t*) is evaluated at the time instants

*t*=

_{k}*k*Δ

*t, k*= 0,1,2,...,

*N*

_{max},

*N*

_{max}

*= t*

_{max}/Δ

*t*, in the time interval [0,

*t*

_{max}]. Since the pulse

*f(t)*is infinitely long in the time domain, it has to be truncated so that only values larger than

*γ*are included into the discrete representation. In many cases, sufficient initial choices are

*γ*= 10

^{−9},

*t*

_{0}= [-ln(

*γ*)]

^{1/2}

*τ*and

*t*

_{max}= 2

*t*

_{0}. Now, the frequency spectrum of

*f*(

*t*) can be calculated via the discrete Fourier transform. The obtained complex frequency spectrum

_{k}*A*(

_{0}*ω*) defines the amplitude and the initial phase for each monochromatic plane wave having angular frequency

_{i}*ω*. Next, the interaction of the monochromatic plane waves with the background structure is solved via the transmittance matrix approach [12]. After the incident fields

_{i}**E**

_{inc}(

*ω*) are solved, the frequency domain source terms

_{i}**S**(

*ω*) are calculated according to Eq. (9). Finally, the time domain source term

_{i}**S**(

*t*) is obtained from the frequency domain terms via the inverse Fourier transform.

When the time domain source term is calculated via the inverse Fourier transform, it is important to verify that spectral sampling resolution, given by *Δf* = 1*/t*_{max}, is adequate. This can be accomplished by calculating the time domain incident electric field **E**_{inc}(*t*) in each layer and checking that the field amplitude has decayed to a negligible value at *t* = *t*_{max}. This guarantees that we have properly sampled the source term in the frequency domain and we obtain an accurate transformation from the frequency to the time domain in Eq. (12). Here, we have assumed that in Eq. (9)
*ε*(*ω*) and *ε*_{inc}(*ω*) are Fourier transformable functions, as is the case with Lorentz, Debye and Drude permittivity distributions.

Probably the most efficient way to evaluate the source term (12) is the fast Fourier transform (FFT). A problem with FFT is that it calculates the source for all time steps *t*_{k}, *k* = 0,1,…,*N*_{max}, at once and thus the computer memory may run out with scatterers consisting of numerous Yee cells. A more viable approach is the discrete Fourier transform that allows the calculation of the source term (12) for each time step at a time. That is,

One way to speed-up the calculation of Eq. (14) is to include only frequencies *ω _{i}* at which the amplitude of the illuminating field |

*A*

_{0}(

*ω*)| is non-negligible.

_{i}## 4. Source term for scatterers in a perfect electric conductor layer

Perfect electric conductor (PEC) is a medium, the electric conductivity of which approaches infinity. Spatially finite PEC scatterers are easily treated in the scattered field technique by applying **E**_{scat} = -**E**_{inc} inside the scatterer [3]. If the stratified background structure of the incident field contains a PEC layer and the PEC layer includes scattering objects extending on its front surface or through the entire PEC layer, a more complicated approach is required to plant the incident field. The source term (12) is not directly evaluable since *ε*_{inc}(*ω*) is not defined for PEC. Thus, we replace the *jωε*_{inc}(*ω*)**E**_{inc}(*ω*) term in Eq. (9) by ∇ × **H**_{inc}(*ω*) according to Eq. (6), and we obtain

Inside a PEC layer, **E**_{inc} = **0** and **H**_{inc} = **0**, which means that also the source term is zero. However, at the surface of the PEC layer fields are non-zero. The boundary conditions [13] state that at the surface of the PEC layer: $\widehat{n}$ × **E**_{inc} = **0**, $\widehat{n}$ × **H**_{inc} = **J**_{s}, and $\widehat{n}$
**⋅ H**_{inc} = 0, where **J**_{s} is the surface current density and $\widehat{n}$ is a unit vector normal to the PEC surface. Since **E**_{inc} ≠ **0** and **H**_{inc} ≠ **0** only on the surface of the PEC layer, the source term is added only to the electric field components that reside on this particular surface. Because tangential components of **H**_{inc} are not continuous across the interface, the ∇ × **H**_{inc}(*ω*) term cannot be evaluated analytically, but we can find a numerical approximation as illustrated next.

Let us consider a PEC interface *z* = *z*_{0} in the *xy*-plane with a unit surface normal $\widehat{n}$ = (0,0,-1). At the interface, *E _{x}*

_{,inc}= 0,

*E*

_{y}_{,inc}= 0 and

*H*

_{z}_{,inc}= 0 due to the boundary conditions. Using the central finite difference approximation, ∇ ×

**H**

_{inc}(

*ω*) term can be expressed at the interface

*z*=

*z*

_{0}as

*z*=

*z*

_{0}+ Δ

*z*/2 are zero since they are located inside the PEC layer. Equation (16.3) can be evaluated analytically. In a typical staggered uniform Yee grid, the electric field components are located in the following coordinate points:

*E*: (

_{x}*i*Δ

*x*, (

*j*+ 1/2)Δ

*y*, (

*k*+ 1/2)Δ

*z*),

*E*: ((

_{y}*i + 1/2)*Δ

*x*,

*j*Δ

*y*, (

*k*+ 1/2)Δ

*z*),

*E*: ((

_{z}*i + 1/2)*Δ

*x*, (

*j*+ 1/2)Δ

*y*,

*k*Δ

*z*), where

*i*,

*j*, and

*k*are integers and Δ

*x,*Δ

*y,*and Δ

*z*define the dimensions of the Yee cell. If we choose that the PEC interface is located at

*z*= (

*k*

_{0}+ 1/2)Δz, then only

*E*and

_{x}*E*components reside on this surface, which means that Eq. (16.3) is not needed at all. Since tangential components of the incident electric field are zero on the PEC interface,

_{y}*x*and

*y*component of the source terms reduces to the form

Since these equations do not contain any frequency dependent multipliers, the time domain representation is obtained by replacing the frequency domain quantities by the corresponding time domain quantities.

## 5. Validation of the scattered field technique

In this section, we verify the validity of the SF technique against the TFSF technique using the two dimensional FDTD method. The FDTD mesh was terminated with the anisotropic perfectly matched absorbing boundaries [14]. Dispersive materials were included into the FDTD model via the auxiliary differential equation [10] and modeled as a single pole Lorentz medium having relative dielectric permittivity of the form:

*ε*

_{s}is the static permittivity,

*ε*

_{∞}is the relative permittivity at infinite frequencies,

*ω*

_{0}is the pole location, and

*δ*

_{0}is the damping factor. In this work, we consider only two dispersive materials, namely silver and gold. Lorentz parameters for these noble metals were found fitting Eq. (18) to the experimental data [15]. For silver, we found the following values provide a good fit:

*ε*= 4.7982x10

_{s}^{7},

*ε*

_{∞}= 4.2167,

*δ*

_{0}= 1.4046×10

^{13}Hz, and

*ω*

_{0}= 2.0254×10

^{12}Hz, and for gold:

*ε*= 2.0397×10

_{s}^{8},

*ε*

_{∞}= 22.3162,

*δ*

_{0}= 4.4183×10

^{13}Hz, and

*ω*

_{0}= 1.2150×10

^{12}Hz.

The SF and TFSF technique are compared by using each to simulate the transmission of a normally incident, TM polarized plane wave pulse through a 100 nm wide slit in a 200 nm thick silver film. The time domain data was transformed to the frequency domain to obtain wavelength dependent transmission efficiency curves. To ensure that the SF technique works properly with dispersive scatterers, the simulations are also performed for the case in which a 80 nm diameter gold cylinder is located in the middle of the slit entrance. In both cases, the slit structure is illuminated by a normally incident plane wave pulse (*τ* = 1.1325 fs, *f*_{0} = *c*_{0}/640nm, *c*_{0} = speed of light in vacuum) propagating in the *z* direction. The TFSF technique was used to model the problem as follows: Since the background structure, where the slit resides, is invariant in the *x* and the *y* directions, propagation of a normally incident plane wave pulse through the background structure can be modeled using a one-dimensional FDTD simulation. This one-dimensional simulation, which is run simultaneously with the actual two-dimensional FDTD simulation, is then used as a look-up table to obtain the incident field for the TFSF technique. The use of look-up tables with the TFSF technique is described in detail in Ref [2].

Light transmission through a slit is usually characterized via the normalized transmission efficiency (*η*), which is the ratio of the transmitted power behind the slit and the incident power at the geometrical entrance of the slit. Figure 1(a)
shows a comparison of the transmission efficiencies obtained by the SF and the TFSF techniques for a 100 nm wide slit in a 200 nm thick silver film as a function of the free space wavelength (λ_{0}). In the case (1), the slit is empty, while in the case (2), a gold cylinder, 80 nm in diameter, is located in the middle of the slit entrance. These results were obtained with a square Yee cell size of 5 nm. Figure 1(b) shows the relative percentage error between the SF and TFSF technique in the case (1) for the square Yee cell sizes of 10, 5, 2.5 and 1 nm. The Fig. 1(c) shows the same data for the case (2). In both cases, good agreement is observed between the techniques. The error clearly diminishes with the decreasing cell size.

## 6. Surface plasmon polariton enhanced light transmission through a slit in a metal film

This section demonstrates the usability of the introduced SF technique in the analysis of surface plasmon polariton (SPP) enhanced light transmittance through a slit in a silver film. A SPP is a composite electromagnetic disturbance that occurs when an incident light field couples resonantly with an oscillating surface charge density distribution. Since the induced surface charge density is proportional to the discontinuity of the electric field component normal to the interface, only transverse magnetic (TM) polarized light can excite SPPs. A semi-infinite dielectric/metal-interface, in which both the metal and the dielectric medium are very thick, supports only a bound SPP mode, the field amplitude of which decays exponentially with distance on the both sides of the interface. SPPs with frequency *ω* have larger momentum than light that propagates in the adjacent dielectric medium with frequency *ω*. Therefore directly incident light cannot excite any SPP modes. Traditionally, SPPs are excited by an evanescent field produced by a prism or a grating that resides close to the dielectric/metal interface.

The inset of Fig. 2(a)
shows a multilayer structure that support excitation of SPPs at the *n*_{2}/*n*_{Ag}-interface. The first layer (*n*_{1} = 2.0) represents the prism layer and the second layer (*n*_{2} = 1.5) is a spacer between the prism and the silver film. The incident TM polarized plane wave having free space wavelength of λ_{0} is incident in the *n*_{1}-layer. The angle of incidence is denoted by *θ*. The refractive index of silver was defined by the single pole Lorentz dispersion model introduced in the previous section. Figure 2(a) shows the calculated reflectance of the structure with *t*_{1} = 480 nm and *t*_{2} = 200 nm. The transmittance through the structure is always zero at the studied wavelength range due to the opaque silver film. When *θ* = 53.7° and λ_{0} = 615 nm, the reflectance of the structure is practically zero and almost all light is coupled to SPPs. The electric field amplitude profiles of *E _{x}* and

*E*components in this resonance case are plotted along the

_{y}*z*-axis in Fig. 2(b). The incident electric field amplitude is 1.0 V/m in the

*n*

_{1}-layer. Due to SPPs, the electric field amplitude is significantly enhanced at the

*n*

_{2}/

*n*

_{Ag}-interface.

Next, we form a 100 nm wide slit into the silver film of the considered multilayer structure and study how SPPs influence on light transmission through the slit. In the SF technique, the first task is to calculate the solution for the incident field in the absence of scatterers. At the SPP resonance angle, a part of the incident plane wave pulse is coupled to SPPs, and a slowly decaying time domain field exists at the *n*_{2}/*n*_{Ag}-interface. This is demonstrated in Fig. 3
, where a green line shows the *H _{y}* component of the illuminating TM-polarized (

*H*,

_{y}*E*,

_{x}*E*) plane wave pulse in the point

_{z}*P*

_{1}at the

*n*

_{1}/

*n*

_{2}-interface. The pulse shape is given by Eq. (13) with

*τ*= 1.6988 fs and

*f*

_{0}=

*c*

_{0}/533.3nm. The red line shows the time domain magnetic field in the point

*P*

_{2}at the surface of the silver film. Clearly, the time duration at which the field decays to zero at the

*n*

_{2}/

*n*

_{Ag}-interface is much larger than the duration of the illuminating pulse due to SPPs propagating along the interface. This has to be taken into account in the spectral sampling resolution (

*Δf*= 1

*/t*

_{max}) of the illuminating plane wave pulse (for details, see Section 3). In practice,

*t*

_{max}has to be so large that at

*t*=

*t*

_{max}field is decayed to zero in all points of the FDTD mesh.

Figure 4
shows the normalized transmission efficiency of a 100 nm wide slit in the silver film of the considered multilayer structure as calculated using the introduced SF technique. The used square Yee cell size was 2.5 nm. With all angles of incidence, the transmitted power through the slit is normalized by the power that would impinge on the slit entrance in the case that *n*_{1} = *n*_{2} = 2.0 and *θ* = 0°. The results clearly show that the excitation of SPPs on the surface of the silver film greatly enhances light transmission through the slit.

## 7. Summary and conclusions

This paper introduces a simple SF technique for dispersive scatterers in dispersive stratified media. An incident time domain plane wave pulse is first decomposed to time harmonic plane waves having angular frequency *ω _{i}*. The interaction of each time harmonic plane wave with the stratified background structure is then solved using the transmittance matrix approach. After the solution for the incident electric fields is found, the frequency domain source term is calculated according to Eq. (9). The same Eq. (9) can be used for all dispersive medium, such as Lorentz, Debye and Drude. Finally, the time domain source term is obtained from the frequency domain terms via the inverse Fourier transform. Dispersive materials are typically incorporated into the FDTD method via the auxiliary differential equation [10] or the recursive convolution method [11]. The introduced technique is compatible with the both approaches. The introduced SF technique is not limited only to plane wave pulses but it can be used with all light beams that can be decomposed to plane waves, as is the case with e.g. the Gaussian and Bessel beams.

The introduced SF technique was verified against the TFSF technique in transmittance of a normally incident plane wave pulse through a 100 nm wide slit in a 200 nm thick silver film. In the second case, a 80 nm diameter gold cylinder was located in the slit entrance. In both cases, a good agreement between the techniques was found. Finally, we studied SPP enhanced light transmittance through a 100nm wide slit in a silver film using the SF technique. The obtained results clearly show that SPPs provide a mechanism to enhance light transmittance through a sub-wavelength slit in a metal film.

## References and links

**1. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media,” IEEE Trans. Antenn. Propag. **14**(3), 302–307 (1966). [CrossRef]

**2. **A. Taflove, and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method* (Second Edition, Artech House, INC., 2000).

**3. **K. S. Kunz, and R. J. Luebbers, *The Finite Difference Time Domain Method for Electromagnetics* (CRC Press, London, 1993).

**4. **S. Winton, P. Kosmas, and C. M. Rappaport, “FDTD simulation of TE and TM plane waves at nonzero incidence in arbitrary layered medium,” IEEE Trans. Antenn. Propag. **53**(5), 1721–1728 (2005). [CrossRef]

**5. **Y.-N. Jiang, D.-B. Ge, and S. J. Ding, “Analysis of TF-SF boundary for 2D-FDTD with plane p-wave propagation in layered dispersive and lossy media,” Prog. In Electromagnetic Res **PIER 83**, 157–172 (2008). [CrossRef]

**6. **J. B. Schneider and K. Abdijalilov, “Analytic field propagation TFSF boundary for FDTD problems involving planar interfaces: PECs, TE, and TM,” IEEE Trans. Antenn. Propag. **54**(9), 2531–2542 (2006). [CrossRef]

**7. **K. Abdijalilov and J. B. Schneider, “Analytic field propagation TFSF boundary for FDTD problems involving planar interfaces: lossy material and evanescent fields,” IEEE Antennas Wirel. Propag. Lett. **5**(1), 454–458 (2006). [CrossRef]

**8. **K. Demarest, R. Plump, and Z. Huan, “FDTD modeling of scatterers in stratified medium,” IEEE Trans. Antenn. Propag. **43**(10), 1164–1168 (1995). [CrossRef]

**9. **S.-C. Kong, J. J. Simpson, and V. Backman, “ADE-FDTD scattered field formulation for dispersive materials,” IEEE Microwave Wireless Comp. Letters **18**(1), 4–6 (2008). [CrossRef]

**10. **M. Okoniewski, M. Mrozowski, and M. A. Stuchly, “Simple treatment of multi-term dispersion in FDTD,” IEEE Microwave Guided Wave Lett. **7**(5), 121–123 (1997). [CrossRef]

**11. **D. F. Kelley and R. J. Luebbers, “Piecewise linear recursive convolution for dispersive media using FDTD,” IEEE Trans. Antenn. Propag. **44**(6), 792–797 (1996). [CrossRef]

**12. **M. G. Moharam, D. A. Pommet, E. B. Grann, and T. K. Gaylord, “Stable implementation of the rigorous coupled-wave analysis for surface relief gratings: enhanced transmittance matrix approach,” J. Opt. Soc. Am. A **12**, 1077–1086 (1995). [CrossRef]

**13. **C.-T. Tai, *Dyadic green functions in electromagnetic theory* (Second Edition, IEEE Press, 1993).

**14. **S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE Trans. Antenn. Propag. **44**(12), 1630–1639 (1996). [CrossRef]

**15. **P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev., B, Solid State **6**(12), 4370–4379 (1972).