Open Access
Add to BibSonomyAbstract
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:
where μ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: where ε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 whereMaxwell’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
whereTraditional total field FDTD kernels designed for dispersive materials can be used to solve discretized versions of Eqs. (10) and (11) with the difference that only the scattered field components are solved instead of the total field components. The term F−1{jωε(ω)Esca(ω)} 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 ∇ × Hsca(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:
where f0 is the central frequency of the pulse, t0 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 tk = kΔt, k = 0,1,2,..., Nmax, Nmax = tmax/Δt, in the time interval [0, tmax]. 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, t0 = [-ln(γ)]1/2 τ and tmax = 2t0. Now, the frequency spectrum of f(tk) can be calculated via the discrete Fourier transform. The obtained complex frequency spectrum A0(ωi) 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 Einc(ωi) are solved, the frequency domain source terms S(ωi) are calculated according to Eq. (9). Finally, the time domain source term 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/tmax, is adequate. This can be accomplished by calculating the time domain incident electric field Einc(t) in each layer and checking that the field amplitude has decayed to a negligible value at t = tmax. 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 tk, k = 0,1,…,Nmax, 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 |A0(ωi)| is non-negligible.
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 Escat = -Einc 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(ω)Einc(ω) term in Eq. (9) by ∇ × Hinc(ω) according to Eq. (6), and we obtain
Inside a PEC layer, Einc = 0 and Hinc = 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: × Einc = 0, × Hinc = Js, and ⋅ Hinc = 0, where Js is the surface current density and is a unit vector normal to the PEC surface. Since Einc ≠ 0 and Hinc ≠ 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 Hinc are not continuous across the interface, the ∇ × Hinc(ω) term cannot be evaluated analytically, but we can find a numerical approximation as illustrated next.
Let us consider a PEC interface z = z0 in the xy-plane with a unit surface normal = (0,0,-1). At the interface, Ex,inc = 0, Ey,inc = 0 and Hz,inc = 0 due to the boundary conditions. Using the central finite difference approximation, ∇ × Hinc(ω) term can be expressed at the interface z = z0 as
where the magnetic field components in the plane z = z0 + Δ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: Ex: (iΔx, (j + 1/2)Δy, (k + 1/2)Δz), Ey: ((i + 1/2)Δx, jΔy, (k + 1/2)Δz), Ez: ((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 = (k0 + 1/2)Δz, then only Ex and Ey 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, x and y component of the source terms reduces to the formSince 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:
where ε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: εs = 4.7982x107, ε∞ = 4.2167, δ0 = 1.4046×1013 Hz, and ω0 = 2.0254×1012 Hz, and for gold: εs = 2.0397×108, ε∞ = 22.3162, δ0 = 4.4183×1013 Hz, and ω0 = 1.2150×1012 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, f0 = c0/640nm, c0 = 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.
Fig. 1 (a) Normalized transmission efficiency of a normally incident plane having free space wavelength of λ0 through a 100 nm wide (w) slit in a 200 nm thick (t) silver film as calculated by the SF and the TFSF technique. In the case (1), the slit is empty, whereas in the case (2), the slit contains a 80 nm diameter gold cylinder in the slit entrance as illustrated in the inset (2). These results were obtained with a square Yee cell size (Δ) of 5 nm. (b) Relative percentage error between the SF and the TFSF techniques as function of the λ0 in the case (1) for square Yee cell sizes of 10, 5, 2.5 and 1 nm. (c) The same as (b), but now a 80 nm diameter gold cylinder is located in the slit entrance.
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 n2/nAg-interface. The first layer (n1 = 2.0) represents the prism layer and the second layer (n2 = 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 n1-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 t1 = 480 nm and t2 = 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 Ex and Ey components in this resonance case are plotted along the z-axis in Fig. 2(b). The incident electric field amplitude is 1.0 V/m in the n1-layer. Due to SPPs, the electric field amplitude is significantly enhanced at the n2/nAg-interface.
Fig. 2 (a) Reflectance of a TM polarized plane wave from the multilayer structure shown in the inset (t1 = 480 nm, t2 = 200 nm, n1 = 2.0, n2 = 1.5, nair = 1.0) as a function of the angle of incidence (θ) and the free space wavelength λ0. The plane wave is incident in the n1-medium. When θ = 53.7° and λ0 = 615 nm, the reflectance of the structure is practically zero and nearly all incident light is coupled to SPPs. At the studied wavelength range, transmittance through the structure is always zero due to the opaque silver film. (b) Amplitude of the electric field components in the studied structure at the SPP resonance angle and the wavelength. Incident electric field amplitude is 1.0 V/m in the n1-medium. The dashed black lines show the locations of the material interfaces of the studied structure.
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 n2/nAg-interface. This is demonstrated in Fig. 3 , where a green line shows the Hy component of the illuminating TM-polarized (Hy, Ex, Ez) plane wave pulse in the point P1 at the n1/n2-interface. The pulse shape is given by Eq. (13) with τ = 1.6988 fs and f0 = c0/533.3nm. The red line shows the time domain magnetic field in the point P2 at the surface of the silver film. Clearly, the time duration at which the field decays to zero at the n2/nAg-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/tmax) of the illuminating plane wave pulse (for details, see Section 3). In practice, tmax has to be so large that at t = tmax field is decayed to zero in all points of the FDTD mesh.
Fig. 3 The structure shown in the inset (n1 = 2.0, n2 = 1.5, t1 = 480nm, t2 = 200nm) is illuminated by a TM polarized (Hy, Ex, Ez) plane wave pulse at the SPP resonance angle (θ = 53.7°). The pulse shape of the illuminating magnetic field is given by Eq. (13) with τ = 1.6988 fs and f0 = c0/533.3nm. The green line shows the magnetic field amplitude of the illuminating field in the point P1, whereas the red line shows Hy in the point P2 at the surface of the silver film. Hy in the point P2 decays slowly to zero due to SPPs propagating along the silver film.
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 n1 = n2 = 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.
Fig. 4 Normalized transmission efficiency of a TM polarized plane wave through a slit structure shown in the inset (n1 = 2.0, n2 = 1.5, t1 = 480nm, t2 = 200nm, w = 100nm) as a function of the free space wavelength λ0. The plane wave is incident in the n1-layer with the incidence angle of θ. When θ = 53.7° and λ0 = 615 nm, light transmission through the slit is enhanced due to the excitation of SPPs.
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).
