OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 17, Iss. 23 — Nov. 9, 2009
  • pp: 21179–21190
« Show journal navigation

Unified perfectly matched layer for finite-difference time-domain modeling of dispersive optical materials

Indika Udagedara, Malin Premaratne, Ivan D. Rukhlenko, Haroldo T. Hattori, and Govind P. Agrawal  »View Author Affiliations


Optics Express, Vol. 17, Issue 23, pp. 21179-21190 (2009)
http://dx.doi.org/10.1364/OE.17.021179


View Full Text Article

Acrobat PDF (272 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Finite-difference time-domain (FDTD) simulations of any electromagnetic problem require truncation of an often-unbounded physical region by an electromagnetically bounded region by deploying an artificial construct known as the perfectly matched layer (PML). As it is not possible to construct a universal PML that is non-reflective for different materials, PMLs that are tailored to a specific problem are required. For example, depending on the number of dispersive materials being truncated at the boundaries of a simulation region, an FDTD code may contain multiple sets of update equations for PML implementations. However, such an approach is prone to introducing coding errors. It also makes it extremely difficult to maintain and upgrade an existing FDTD code. In this paper, we solve this problem by developing a new, unified PML algorithm that can effectively truncate all types of linearly dispersive materials. The unification of the algorithm is achieved by employing a general form of the medium permittivity that includes three types of dielectric response functions, known as the Debye, Lorentz, and Drude response functions, as particular cases. We demonstrate the versatility and flexibility of the new formulation by implementing a single FDTD code to simulate absorption of electromagnetic pulse inside a medium that is adjacent to dispersive materials described by different dispersion models. The proposed algorithm can also be used for simulations of optical phenomena in metamaterials and materials exhibiting negative refractive indices.

© 2009 Optical Society of America

1. Introduction

Contemporary optical problems often require full-blown, three-dimensional numerical simulations carried out using the finite-difference time-domain (FDTD) technique [1

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

, 2

2. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, 2005).

]. Since all real materials exhibit dispersion (frequency dependence) of their optical properties, physical models adequately representing their response to electromagnetic fields have to be incorporated into FDTD simulators. In any FDTD code, one often needs to make an important tradeoff between numerical stability and computational cost [3

3. P. Petropoulos, “Stability and phase error analysis of FD-TD in dispersive dielectrics,” IEEE Trans. Antennas Propag. 42, 62–69 ( 1994). [CrossRef]

, 4

4. M. Premaratne and S. K. Halgamuge, “Rigorous analysis of numerical phase and group velocity bounds in Yee’s FDTD grid,” IEEE Microwave Wirel. Compon. Lett. 17, 556–558 ( 2007). [CrossRef]

]. Consequently, different compromises result in different variants of FDTD formulations for a particular model of material dispersion. Among the widely used variants are: the recursive convolution (RC) method [5

5. R. Luebbers, F. Hunsberger, K. Kunz, R. Standler, and M. Schneider, “A frequency-dependent finite-difference time-domain formulation for dispersive materials,” IEEE Trans. Electromagn. Compat. 32, 222–227 ( 1990). [CrossRef]

], the piecewise linear recursive convolution (PLRC) method [6

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

], the auxiliary differential equation (ADE) method [7

7. J. Young, “Propagation in linear dispersive media: Finite difference time-domain methodologies,” IEEE Trans. Antennas Propag. 43, 422–426 ( 1995). [CrossRef]

], and the Z-transform method [8

8. D. Sullivan, “Z-transform theory and the FDTD method,” IEEE Trans. Antennas Propag. 44, 28–34 ( 1996). [CrossRef]

]. Unfortunately, none of these methods is flexible enough to provide optimum accuracy for arbitrarily dispersive materials, and they demand a tailor-made implementation for each specific case. For example, in the ADE method, FDTD formulations for the Debye and Lorentz models of dielectric response are different, thus requiring re-derivation of update equations whenever materials with complex dispersion models need to be simulated [9

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

, 10

10. Y. Takayama and W. Klaus, “Reinterpretation of the auxiliary differential equation method for FDTD,” IEEE Microwave Wirel. Compon. Lett. 12, 102–104 ( 2002). [CrossRef]

]. Obviously, such formulations are inefficient since they may inadvertently introduce errors into an existing FDTD code.

Recently, Han et al. showed that dispersive materials described by the Debye and Lorentz models can be represented by a generic expression template consisting of complex-conjugate pairs of simple residues [11

11. M. Han, R. Dutton, and S. Fan, “Model dispersive media in finite-difference time-domain method with complex-conjugate pole-residue pairs,” IEEE Microwave Compon. Lett. 16, 119–121 ( 2006). [CrossRef]

]. Adoption of this novel representation enables one to create a single FDTD code for all practically useful media and eliminate the deficiencies associated with multiple formulations. However, most of these simulations require truncation of an unbounded FDTD simulation region by an electromagnetically bounded region, which can only be achieved by having perfectly-matched layers (PMLs) tailor-made for the dispersive materials of interest. Therefore, depending on the number of dispersive materials being truncated at boundaries of a simulation region, the corresponding FDTD code may contain multiple sets of update equations for implementing PMLs. Such a tailor-made implementation, and its sub-sequent integration with the existing FDTD code, is prone to introduce coding errors and also makes the maintenance and upgrading of the code extremely difficult. In this paper, for the first time to the best of our knowledge, we derive an algorithm for a unified PML and demonstrate its practical importance through a number of numerical simulations. In contrast to previously reported algorithms, based on different auxiliary equations for each type of media [12

12. K. P. Prokopidis, “On the development of efficient FDTD-PML formulations for general dispersive media,” Int. J. Numer. Model. 21, 395–411 ( 2008). [CrossRef]

], the new approach allows us to create a single FDTD code that can be used for a variety of optical simulations requiring suppression of reflections from the grid boundaries.

This paper is organized as follows. We start in Section 2 by introducing the generalized expression for relative permittivity of an arbitrary linear optical media. Using the methodology of complex-frequency-shifted PML, we then derive the relevant update equations for auxiliary functions which are independent of the PML parameters. These parameters are included in the FDTD implementation by restoring the electromagnetic field from auxiliary functions in the PML region. The restoration procedure and the steps required to update the values of electro-magnetic field on the FDTD grid are given at the end of Section 2. In Section 3, we illustrate the flexibility of an FDTD code based on the proposed algorithm by applying it to a number of realistic situations. We summarize our results and conclude the paper in Section 4.

2. Construction of unified PML for arbitrary dispersive media

To develop a unified PML capable of truncating different dispersive media, we need to represent the dielectric permittivities of these media in a generalized form. Since the optical media are nonmagnetic in most situations, we neglect dispersion of the magnetic permeability, µ, throughout this study. According to Han et al. [11

11. M. Han, R. Dutton, and S. Fan, “Model dispersive media in finite-difference time-domain method with complex-conjugate pole-residue pairs,” IEEE Microwave Compon. Lett. 16, 119–121 ( 2006). [CrossRef]

], the frequency dependence of relative permittivity of any dielectric medium can be expanded using a set of complex parameters ap and cp as follows:

ε˜(ω)=ε+p=0p(cpiωap+cp*iωap*),
(1)

Table 1. Poles and their residues for CCPR media corresponding to Debye, Lorentz, and Drude models of the dielectric response

table-icon
View This Table
| View All Tables

aΔεp is the weight of the pth pole and τp is the carrier relaxation time.

b ωp is the resonance frequency and δp is the damping parameter; parameters a 0 and c 0 are the same as for Debye model.

c ω pl and ωc are the plasma and collision frequencies, respectively.

where ε is the real-valued, high-frequency permittivity, P is a positive integer, and an asterisk denotes complex conjugation. Hereafter, we use a tilde above a time-domain function to denote its Fourier transform in the frequency domain, i.e.,

A˜(ω)=+A(t)exp(iωt)dt.

The values of poles, ap, and residues, cp, entering Eq. (1) depend on a particular model of dispersive media. For generality of this treatment, we refer to materials that are described by Eq. (1) as complex-conjugate-pole residue (CCPR) medium from here onwards. Table 1 provides expressions of permittivities and summarizes values of the parameters ap and cp for three widely used models of dielectric response from dispersive media [2

2. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, 2005).

, 13

13. R. Luebbers, F. Hunsberger, and K. Kunz, “A frequency-dependent finite-difference time-domain formulation for transient propagation in plasma,” IEEE Trans. Antennas Propag. 39, 29–34 ( 1991). [CrossRef]

]. In the case of Debye and Lorentz models, we have explicitly included the static electrical conductivity σ0=- 0 limω→0 ωε̃ (ω) into the expressions of dielectric permittivities. The Drude model does not require this modification as it accounts for σ00 ω 2 pl ω -1 c automatically.

2.1. Synthesis of PML for CCPR media

Let Ẽ(r,ω) and H̃(r,ω) be the electric and magnetic fields inside a CCPR medium. In the absence of external sources, the evolution of these fields (in the frequency domain) is governed by the Maxwell equations,

×E˜=iωμ0μH˜,
(2a)
×H˜=iωε0ε˜(ω)E˜,
(2b)

where µ 0 and ε 0 are the magnetic and dielectric permittivities of free space. An absorbing PML that is perfectly impedance-matched to its adjacent media, can be created by introducing a specific transformation of the Cartesian coordinates (for details see [14

14. W. C. Chew and W. H. Weedon, “A 3-D perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microwave Opt. Tech. Lett 7, 599–604 ( 1994). [CrossRef]

, 15

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

]). With this transformation, in the PML region Eqs. (2) take the form

'×E˜=iωμ0μΓ˜H˜,
(3a)
'×H˜=iωε0ε˜(ω)Γ˜E˜,
(3b)

where ∇′ is the transformed operator and the PML material tensor is given by Γ˜=(syszsx000sxszsy000sxsysz).

The presence of material tensor Γ˜ on the right-hand side of Eqs. (3) is equivalent to filling of PML by an artificial absorbing material with anisotropic permittivity ε̃ (ω) Γ˜ and anisotropic permeability µ Γ˜ . Clearly, for any choice of the stretching coefficients sj (j=x,y,z), the impedances of CCPR medium and PML are equal with the value η=μ0μ(ε0ε˜). Attenuation properties of the PML in the directions of the spatial coordinates can be controlled by choosing the coefficients s j appropriately. In the present paper, we use the following definition [16

16. M. Kuzuoglu and R. Mittra, “Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers,” IEEE Microwave Guided Wave Lett. 6, 447–449 ( 1996). [CrossRef]

, 17

17. S. Gedney, “Scaled CFS-PML: It is more robust, more accurate, more efficient, and simple to implement. Why aren’t you using it?” in Antennas and Propagation Society International Symposium364–367 ( 2005).

]:

sj=κj+σjγiωε0.
(4)

Such a PML is referred to as a complex-frequency-shifted (CFS) PML and is renowned for its improved absorption efficiency as compared to other PML models [19

19. J.-P. Berenger, “Numerical reflection from FDTD-PMLs: A comparison of the split PML with the unsplit and CFS PMLs,” IEEE Trans. Antennas Propag. 50, 258–265 ( 2002). [CrossRef]

, 18

18. D. Correia and J.-M. Jin, “Performance of regular PML, CFS-PML, and second-order PML for waveguide problems,” Microwave Opt. Technol. Lett. 48, 2121–2126 ( 2006). [CrossRef]

]. The impact of parameters σj>0, κj>1, and γ>0 on the attenuation properties of CFS PML is thoroughly analyzed in Refs. [19

19. J.-P. Berenger, “Numerical reflection from FDTD-PMLs: A comparison of the split PML with the unsplit and CFS PMLs,” IEEE Trans. Antennas Propag. 50, 258–265 ( 2002). [CrossRef]

, 21

21. J.-P. Berenger, “Application of the CFS PML to the absorption of evanescent waves in waveguides,” IEEE Microwave Wirel. Compon. Lett. 12, 218–220 ( 2002). [CrossRef]

, 20

20. D. Correia and J.-M. Jin, “On the development of a higher-order PML,” IEEE Trans. Antennas Propag. 53, 4157–4163 ( 2005). [CrossRef]

]. From a physical point of view, σj resembles the conductivity profile of PML and provides attenuation of propagating waves in the jth direction, whereas κj and γ absorb evanescent waves. To avoid numerical reflections from the PML boundaries, the parameters σj and κj need to be smoothly varied in space. One choice of these parameters and the guidelines for the coefficients involved can be found in [22

22. Y. Rickard and N. Georgieva, “Problem-independent enhancement of PML ABC for the FDTD method,” IEEE Trans. Antennas Propag. 51, 3002–3006( 2003). [CrossRef]

]:

σj=σmax(k/δj)m+n,           κj=1+(κmax1)(k/δj)n,

where kδj, δj is the PML depth in the jth direction, m ∈ [-3,3] and n ∈ [2,6] are the user-defined integers meeting the condition m+n > 1, κ max ∈ (1,10], σmax=- 0(m+n+1) lnR 0/(2δj), c is the speed of light in vacuum, and R 0 ∈ [10-12,10-2].

2.2. Derivation of FDTD update equations

Implementation of Eqs. (2) and (3) using the FDTD algorithm requires discretization of the electric and magnetic fields on the FDTD grid and derivation of the update equations for them. It is desirable that update equations do not depend on the exact form of PML material tensor because such a dependance requires their modification for each specific choice of the attenuation parameters sj. This can be realized by rewriting the update equations in terms of two auxiliary functions,

R˜E=Γ˜E˜,R˜H=Γ˜H˜.
(5)

Once values of these functions are found in the time domain, they are used to calculate the real fields, E(r, t) and H(r, t).

The update equations for R E and R H can be derived as follows. Substituting Eq. (1) into Eq. (3b) and omitting the dash sign for convenience, we are led to the equation

×H˜=iωε0εR˜E+p=0P(J˜p+K˜p),
(6)

where

J˜p=iωε0cpiωapR˜E,K˜p=iωε0cp*iωap*R˜E.
(7)

Discretizing Eqs. (7) in the time domain on the FDTD grid, we obtain

Jpn+1=αpJpn+βp(REn+1REn),Kpn=(Jpn)*,
(8)

where αp=1apΔt21+apΔt2,βp=ε0cp1+apΔt2.

Using Eqs. (8), we can write the discrete time-domain version of Eq. (6) in the form ×Hn=ε0εREn+12REn12Δt +p=0P[(1+αp)Jpn12]+(REn+12REn12)p=0P(βp),

where ℜ sign stands for real part of a complex number. From this equation, we find the following update equation for R E:

REn+12=REn12+×Hn+[(1+αp)Jpn12]ε0εΔt(βp).
(9)

Similarly, using Eq. (3a), it is easy to show that the update equation for RH is:

RHn+1=RHn×En+12μ0μΔt.
(10)

2.3. Restoration of the electromagnetic field inside PML

In order to restore the electric and magnetic fields from auxiliary functions, we need to find the relationships between them in discrete time-domain. This can be done by converting Eqs. (5) to continuous time-domain using the Fourier operator pair -/∂t and discretizing the resulting differential equations with appropriate difference operators. This procedure, however, is quite tedious in our case and becomes even harder for more complicated forms of the coefficients sj (e.g., for negative-index-materials [23

23. S. Cummer, “Perfectly matched layer behavior in negative refractive index materials,” IEEE Antennas Wirel. Propag. Lett. 3, 172–175 ( 2004). [CrossRef]

] or higher-order PMLs [20

20. D. Correia and J.-M. Jin, “On the development of a higher-order PML,” IEEE Trans. Antennas Propag. 53, 4157–4163 ( 2005). [CrossRef]

]). The inconvenience associated with discretizing of a high-order differential equation can be avoided by employing the Z-transform technique [24

24. A. V. Oppenheim and R. W. Schafer, Digital Signal Processing (Prentice-Hall, 1975).

, 8

8. D. Sullivan, “Z-transform theory and the FDTD method,” IEEE Trans. Antennas Propag. 44, 28–34 ( 1996). [CrossRef]

]. We use this technique in the following formulation.

Without loss of generality, let us derive the difference equation for the x component of the electric field using the definition, R̃E,x=Γ˜ x Ẽx. Introducing new parameters, ξj=ξ 0j/(κj ε 0), ξ 0=γ/ε 0, and utilizing Eq. (4), the relation between R̃E,x and Ẽx takes the form κxκyκzR˜E,x(ξy)(ξz)=E˜x(ξx)(ξ0).

Recasting this equation into the Z-domain using the transform pair 1ξ11z1eξΔt,

we find κxκyκzRE,x(z)1z1(eξyΔt+eξzΔt)+z2e(ξy+ξz)Δt=Ex(z)1z1(eξxΔt+eξ0Δt)+z2e(ξx+ξ0)Δt.

This equation can be readily inverted to obtain its equivalent time-domain form suitable for numerical implementation. Using the transform pair z -k An(z)↔A n-k, we get

Exn+1=Exn(eξyΔt+eξzΔt)Exn1e(ξy+ξz)Δt
+κxκyκz[RE,xn+1RE,xn(eξxΔt+eξ0Δt)+RE,xn1e(ξx+ξ0)Δt].
(11)

Analogous expressions can be obtained for the Ey and Ez components from Eq. (11) upon index permutations xy and xz, respectively. The components H n+1j satisfy equations identical to those for E n+1 j if we use Rk H,j in place of Rk E,j.

2.4. Steps to update FDTD grid

The equations derived in the preceding subsection allow us to come up with an unified FDTD algorithm for simulation of optical phenomena in a wide range of dispersive media truncated by a perfectly matched absorbing boundary. The main steps of the algorithm for updating the electromagnetic field in PML can be given as

1. Apply the update equation (9) and store R n+1/2 E values. The back-stored values of H n, J n-1/2 p, and R n-1/2 E are required for this update.

2. Calculate and store the auxiliary variables J n+1/2 p for each pole p using Eq. (8). The value of R n+1/2 E are then calculated for the current time step using the back-stored values of J n-1/2 p and R n-1/2 E.

3. Use Eq. (11) and similar equations for y and z components to restore the electric field E n+1/2 from the back-stored values of E n-1/2, E n-3/2, R n-1/2 E, and R n-3/2 E as well as from R n+1/2 E calculated during the current time step.

4. Update R n+1 H values using Eq. (10). The back-stored values of R nH and E n+1/2 are required for this update.

5. Employ an analog of Eq. (11) to calculate H n+1 using back-stored values of H n, H n-1, R nH, and R n-1 H as well as the value of R n+1 H calculated during the current time step.

3. Numerical examples

In this section, we illustrate usefulness of the developed unified PML for several scenarios of TM-wave propagation in a two-dimensional FDTD grid. As a measure of the effectiveness

Table 2. Material parameters used in the simulations

table-icon
View This Table
| View All Tables

of the PML, we use deviation of the electric field component Ez, calculated using a 50×50-cell grid truncated with a PML, from that calculated on a reference grid with 400×400 cells. Mathematically, this deviation may be characterized by the global error [12

12. K. P. Prokopidis, “On the development of efficient FDTD-PML formulations for general dispersive media,” Int. J. Numer. Model. 21, 395–411 ( 2008). [CrossRef]

, 25

25. J.-P. Berenger, “Three-dimensional perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 127, 363–379 ( 1996). [CrossRef]

]

χ2=x=150y=150[Ez(x,y)Eref,z(x,y)]2.

The duration of simulations is chosen to ensure that electromagnetic waves do not experience reflections from the boundaries of the reference grid. To excite the TM wave, we place a hard electric source at the center of the simulation region and produce a Gaussian pulse with the electric field

Ez(n)=exp[(n5010)2]sin(2πfcnΔt),
(12)

where n is the time-step number and fc is the carrier frequency.

The thickness of the PML is assumed to be equal in the x and y directions and its parameters (see page 5) are empirically chosen to be m=0, n=4, κ max=2, γ=1, and R 0=10-7. The cell dimensions Δxy=λ min/40, where λ min=c/(2.5 fc) is the wavelength of interest, and the Courant factor S=cΔtx=0.1.

We examine the PML efficiency using several types of dispersive materials representing different models of dielectric response. More specifically, we consider a 6-pole Debye model of human muscles [26

26. O. P. Gandhi, B.-Q. Gao, and J.-Y. Chen, “A frequency-dependent finite-difference time-domain formulation for general dispersive media,” IEEE Trans. Microwave Theory Tech. 41, 658–665 ( 1993). [CrossRef]

, 27

27. W. D. Hurt, “Multiterm Debye dispersion relations for permittivity of muscle,” IEEE Trans. Bio. Eng. 32, 60–64 ( 1985). [CrossRef]

], a 2-pole Lorentz model of fluoride glass [28

28. K. E. Oughstun, “Computational methods in ultrafast time-domain optics,” Comput. Sci. Eng. 5, 22–32 ( 2003). [CrossRef]

], and a 2-pole Drude model of cold plasma [12

12. K. P. Prokopidis, “On the development of efficient FDTD-PML formulations for general dispersive media,” Int. J. Numer. Model. 21, 395–411 ( 2008). [CrossRef]

]. It is also interesting to look at the material whose dielectric response cannot be adequately described by these models and requires a more complex representation. As an example of such medium, we consider gold. The Drude model of gold is known to be unsuitable for energies greater than 1.8 eV. A more accurate dispersion of gold’s dielectric permittivity is given by the Drude–Lorentz equation [29

29. J. Xi and M. Premaratne, “Analysis of the optical force dependency on beam polarization: Dielectric/metallic spherical shell in a Gaussian beam,” J. Opt. Soc. Am. B 26, 973–980 ( 2009). [CrossRef]

, 30

30. A. Vial, A.-S. Grimault, D. Macías, D. Barchiesi, and M. L. 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(1–7) ( 2005). [CrossRef]

]

ε˜(ω)=εωpl2ω(ω+iωc)+Δε1ω12ω122iωδ1ω2.
(13)

The material parameters characterizing the above four media are summarized in Table 2, and the corresponding permittivity spectra are illustrated by Fig. 1.

Fig. 1. Frequency dependence of the dielectric permittivity of four different materials using (a) a 6-pole Debye, (b) a 2-pole Lorentz, (c) a 2-pole Drude, and (d) a 3-pole Drude–Lorentz models. See Table 2 for material parameters.

The FDTD algorithm of Section 2 allows us to simulate electromagnetic pulse propagation in all four media by employing a single FDTD code. To utilize this code, we calculate the values of ap and cp using Tables 1 and 2. When only the experimental spectrum of dielectric response function is known, one needs to employ curve fitting using Eq. (1). In the simulations, we chose the carrier frequency of the excitation pulse to be 1 GHz for the Debye and Drude media and 500 THz for the Lorentz and Drude–Lorentz media. Figure 2 shows evolution of global error with time for sample media truncated by PMLs of different thicknesses. As seen there, the global error is relatively small in the case of an 8-cell PML (blue curves) for all four dispersive media. In particular, the maximum value of χ 2 is at most ~10-3, and can be ~10-6 for the Drude model, signifying high efficiency of the PML (for comparison, see Ref. [12

12. K. P. Prokopidis, “On the development of efficient FDTD-PML formulations for general dispersive media,” Int. J. Numer. Model. 21, 395–411 ( 2008). [CrossRef]

]). Obviously, thinner PMLs (red and green curves) result in bigger global errors since the waves are attenuated less if they travel a smaller distance within the PML. It should be also noted that the time required for saturation of global error depends on the group velocity at the pulse carrier frequency and, therefore, different for all four media.

In practice, one often needs to model propagation of optical waves through complicated heterostructures containing different types of dispersive materials. If dimensions of the heterostructures exceed a computationally manageable grid, PMLs to each of the bulky materials need to be employed to limit the simulation region. As explained earlier, the proposed algorithm allows one to naturally include any number of different PMLs in one FDTD code by only setting a proper map of material parameters. To illustrate this possibility, we consider penetration of light from vacuum into two half-spaces filled up with fluoride glass on left and with gold on right (for parameters, see Table 2). As before, an excitation pulse in the form of Eq. (12) with fc=500 THz is launched at the center of 380×380-cell grid surrounded by a 10-cell PML.

Figure 3 shows the snapshots of electric field in the whole grid at the time steps n=2300 (left panel) and n=3300 (right panel). One can see that reflections from all parts of the PML are almost negligible for both the oblique and grazing waves. It should also be noted that the PML provides efficient absorption for both the high-frequency (left panel) and low-frequency (right panel) field components. These conclusions are supported by the evolution of global error shown by the red curve in Fig. 4. Particularly, χ 2≈3×10-8 for the left panel in Fig. 3, but it becomes ≈1.2×10-3 for the right panel in Fig. 3. As expected, the global error drastically depends on the PML thickness, increasing by more than a factor of 20 for a 5-cell PML (green curve in Fig. 4) and decreasing nearly by a factor of 6 for a 15-cell PML (blue curve in Fig. 4) at n=3300.

Fig. 2. Evolution of global error for four dispersive media shown in Fig. 1: (a) 5-pole Debye medium; (b) 3-pole Lorentz medium; (c) 2-pole Drude medium; (d) 3-pole Drude–Lorentz medium. The inset in panel (a) shows schematically three different PMLs surrounding the simulation region. For simulation parameters, see the text.

Even though preceding simulations were restricted to two-dimensional waves excited by point sources, they can be easily repeated in a three-dimensional grid with other types of sources, e.g., with a plane-wave source realized by the total-field-scattered-field technique [2

2. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, 2005).

]. These and other numerical examples indicate that the proposed algorithm of CFS PML provides an efficient absorbing boundary for CCPR media, and therefore it is suitable for any dispersive material exhibiting linear response to an optical field.

4. Conclusion

Owing to the simplicity and stability of the standard FDTD algorithm, it is widely used for simulating complex optical structures made of different dispersive materials. To reduce the computational time, it is common to limit the size of the FDTD grid using artificial, perfectly matched layers. In this paper, we highlighted some of the existing strategies for devising PMLs for truncating FDTD grids. However, each of these implementations add an added layer of complexity that makes the coding and maintenance of the FDTD algorithm error prone and complex. To address this issue, we presented in this paper a generalized algorithm for implementing a unified PML model that can be applied to a variety of dispersive media widely used in practice. Based on this algorithm, we developed an optimized FDTD code suitable for simulations of infinite and semi-infinite heterostructures containing any number of dispersive materials described by the Debye, Lorentz, and Drude models, or any other model that could be expressed as a linear combination of these models (e.g., gold). We demonstrated the effectiveness of the unified PML by applying it to several problems involving multiple materials with different dispersion properties. We believe that the proposed algorithm will be beneficial for researches who use the FDTD scheme for a direct solution of the electromagnetic field equations.

Fig. 3. Density plots of the electric field component, Ez, after 2300 (left panel) and 3300 (right panel) FDTD time-steps. A Gaussian pulse, emitted by a point source (cross) located in vacuum (region C) illuminates region A containing fluoride glass and the region B containing gold. The parts of a 10-cell PML marked as A′, B′, and C′ are matched to the regions A, B, and C, respectively. For simulation parameters, see the text.
Fig. 4. Evolution of global error with the time-step number during simulations shown in Fig. 3 (red curve). For comparison, blue and green curves show error evolution for thicker and thinner PMLs. All three curves are calculated for a 380×380 grid and compared to a 1000×1000 reference grid.

Acknowledgments

This work was funded by the Australian Research Council through its Discovery Grant scheme under grant DP0877232. The work of GPA is supported by the NSF award ECCS-0801772.

References and links

1.

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

2.

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, 2005).

3.

P. Petropoulos, “Stability and phase error analysis of FD-TD in dispersive dielectrics,” IEEE Trans. Antennas Propag. 42, 62–69 ( 1994). [CrossRef]

4.

M. Premaratne and S. K. Halgamuge, “Rigorous analysis of numerical phase and group velocity bounds in Yee’s FDTD grid,” IEEE Microwave Wirel. Compon. Lett. 17, 556–558 ( 2007). [CrossRef]

5.

R. Luebbers, F. Hunsberger, K. Kunz, R. Standler, and M. Schneider, “A frequency-dependent finite-difference time-domain formulation for dispersive materials,” IEEE Trans. Electromagn. Compat. 32, 222–227 ( 1990). [CrossRef]

6.

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

7.

J. Young, “Propagation in linear dispersive media: Finite difference time-domain methodologies,” IEEE Trans. Antennas Propag. 43, 422–426 ( 1995). [CrossRef]

8.

D. Sullivan, “Z-transform theory and the FDTD method,” IEEE Trans. Antennas Propag. 44, 28–34 ( 1996). [CrossRef]

9.

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

10.

Y. Takayama and W. Klaus, “Reinterpretation of the auxiliary differential equation method for FDTD,” IEEE Microwave Wirel. Compon. Lett. 12, 102–104 ( 2002). [CrossRef]

11.

M. Han, R. Dutton, and S. Fan, “Model dispersive media in finite-difference time-domain method with complex-conjugate pole-residue pairs,” IEEE Microwave Compon. Lett. 16, 119–121 ( 2006). [CrossRef]

12.

K. P. Prokopidis, “On the development of efficient FDTD-PML formulations for general dispersive media,” Int. J. Numer. Model. 21, 395–411 ( 2008). [CrossRef]

13.

R. Luebbers, F. Hunsberger, and K. Kunz, “A frequency-dependent finite-difference time-domain formulation for transient propagation in plasma,” IEEE Trans. Antennas Propag. 39, 29–34 ( 1991). [CrossRef]

14.

W. C. Chew and W. H. Weedon, “A 3-D perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microwave Opt. Tech. Lett 7, 599–604 ( 1994). [CrossRef]

15.

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

16.

M. Kuzuoglu and R. Mittra, “Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers,” IEEE Microwave Guided Wave Lett. 6, 447–449 ( 1996). [CrossRef]

17.

S. Gedney, “Scaled CFS-PML: It is more robust, more accurate, more efficient, and simple to implement. Why aren’t you using it?” in Antennas and Propagation Society International Symposium364–367 ( 2005).

18.

D. Correia and J.-M. Jin, “Performance of regular PML, CFS-PML, and second-order PML for waveguide problems,” Microwave Opt. Technol. Lett. 48, 2121–2126 ( 2006). [CrossRef]

19.

J.-P. Berenger, “Numerical reflection from FDTD-PMLs: A comparison of the split PML with the unsplit and CFS PMLs,” IEEE Trans. Antennas Propag. 50, 258–265 ( 2002). [CrossRef]

20.

D. Correia and J.-M. Jin, “On the development of a higher-order PML,” IEEE Trans. Antennas Propag. 53, 4157–4163 ( 2005). [CrossRef]

21.

J.-P. Berenger, “Application of the CFS PML to the absorption of evanescent waves in waveguides,” IEEE Microwave Wirel. Compon. Lett. 12, 218–220 ( 2002). [CrossRef]

22.

Y. Rickard and N. Georgieva, “Problem-independent enhancement of PML ABC for the FDTD method,” IEEE Trans. Antennas Propag. 51, 3002–3006( 2003). [CrossRef]

23.

S. Cummer, “Perfectly matched layer behavior in negative refractive index materials,” IEEE Antennas Wirel. Propag. Lett. 3, 172–175 ( 2004). [CrossRef]

24.

A. V. Oppenheim and R. W. Schafer, Digital Signal Processing (Prentice-Hall, 1975).

25.

J.-P. Berenger, “Three-dimensional perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 127, 363–379 ( 1996). [CrossRef]

26.

O. P. Gandhi, B.-Q. Gao, and J.-Y. Chen, “A frequency-dependent finite-difference time-domain formulation for general dispersive media,” IEEE Trans. Microwave Theory Tech. 41, 658–665 ( 1993). [CrossRef]

27.

W. D. Hurt, “Multiterm Debye dispersion relations for permittivity of muscle,” IEEE Trans. Bio. Eng. 32, 60–64 ( 1985). [CrossRef]

28.

K. E. Oughstun, “Computational methods in ultrafast time-domain optics,” Comput. Sci. Eng. 5, 22–32 ( 2003). [CrossRef]

29.

J. Xi and M. Premaratne, “Analysis of the optical force dependency on beam polarization: Dielectric/metallic spherical shell in a Gaussian beam,” J. Opt. Soc. Am. B 26, 973–980 ( 2009). [CrossRef]

30.

A. Vial, A.-S. Grimault, D. Macías, D. Barchiesi, and M. L. 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(1–7) ( 2005). [CrossRef]

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(160.4670) Materials : Optical materials
(260.2030) Physical optics : Dispersion
(080.1753) Geometric optics : Computation methods
(050.1755) Diffraction and gratings : Computational electromagnetic methods

ToC Category:
Physical Optics

History
Original Manuscript: September 23, 2009
Revised Manuscript: October 23, 2009
Manuscript Accepted: November 4, 2009
Published: November 6, 2009

Citation
Indika Udagedara, Malin Premaratne, Ivan D. Rukhlenko, Haroldo T. Hattori, and Govind P. Agrawal, "Unified perfectly matched layer for finite-difference time-domain modeling of dispersive optical materials," Opt. Express 17, 21179-21190 (2009)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-23-21179


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media," IEEE Trans. Antennas Propag. 14, 302-307 (1966). [CrossRef]
  2. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, 2005).
  3. P. Petropoulos, "Stability and phase error analysis of FD-TD in dispersive dielectrics," IEEE Trans. Antennas Propag. 42, 62-69 (1994). [CrossRef]
  4. M. Premaratne and S. K. Halgamuge, "Rigorous analysis of numerical phase and group velocity bounds in Yee’s FDTD grid," IEEE Microwave Wirel. Compon. Lett. 17, 556-558 (2007). [CrossRef]
  5. R. Luebbers, F. Hunsberger, K. Kunz, R. Standler, and M. Schneider, "A frequency-dependent finite-difference time-domain formulation for dispersive materials," IEEE Trans. Electromagn. Compat. 32, 222-227 (1990). [CrossRef]
  6. D. Kelley and R. Luebbers, "Piecewise linear recursive convolution for dispersive media using FDTD," IEEE Trans. Antennas Propag. 44, 792-797 (1996). [CrossRef]
  7. J. Young, "Propagation in linear dispersive media: Finite difference time-domain methodologies," IEEE Trans. Antennas Propag. 43, 422-426 (1995). [CrossRef]
  8. D. Sullivan, "Z-transform theory and the FDTD method," IEEE Trans. Antennas Propag. 44, 28-34 (1996). [CrossRef]
  9. M. Okoniewski, M. Mrozowski, and M. Stuchly, "Simple treatment of multi-term dispersion in FDTD," IEEE Microwave Guided Wave Lett. 7, 121-123 (1997). [CrossRef]
  10. Y. Takayama and W. Klaus, "Reinterpretation of the auxiliary differential equation method for FDTD," IEEE Microwave Wirel. Compon. Lett. 12, 102-104 (2002). [CrossRef]
  11. M. Han, R. Dutton, and S. Fan, "Model dispersive media in finite-difference time-domain method with complexconjugate pole-residue pairs," IEEE Microwave Compon. Lett. 16, 119-121 (2006). [CrossRef]
  12. K. P. Prokopidis, "On the development of efficient FDTD-PML formulations for general dispersive media," Int. J. Numer. Model. 21, 395-411 (2008). [CrossRef]
  13. R. Luebbers, F. Hunsberger, and K. Kunz, "A frequency-dependent finite-difference time-domain formulation for transient propagation in plasma," IEEE Trans. Antennas Propag. 39, 29-34 (1991). [CrossRef]
  14. W. C. Chew and W. H. Weedon, "A 3-D perfectly matched medium from modified Maxwell’s equations with stretched coordinates," Microwave Opt. Tech. Lett 7, 599-604 (1994). [CrossRef]
  15. S. Gedney, "An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices," IEEE Trans. Antennas Propag. 44, 1630-1639 (1996). [CrossRef]
  16. M. Kuzuoglu and R. Mittra, "Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers," IEEE Microwave Guided Wave Lett. 6, 447-449 (1996). [CrossRef]
  17. S. Gedney, "Scaled CFS-PML: It is more robust, more accurate, more efficient, and simple to implement. Why aren’t you using it?" in Antennas and Propagation Society International Symposium 364-367 (2005).
  18. D. Correia and J.-M. Jin, "Performance of regular PML, CFS-PML, and second-order PML for waveguide problems," Microwave Opt. Technol. Lett. 48, 2121-2126 (2006). [CrossRef]
  19. J.-P. Berenger, "Numerical reflection from FDTD-PMLs: A comparison of the split PML with the unsplit and CFS PMLs," IEEE Trans. Antennas Propag. 50, 258-265 (2002). [CrossRef]
  20. D. Correia and J.-M. Jin, "On the development of a higher-order PML," IEEE Trans. Antennas Propag. 53, 4157-4163 (2005). [CrossRef]
  21. J.-P. Berenger, "Application of the CFS PML to the absorption of evanescent waves in waveguides," IEEE Microwave Wirel. Compon. Lett. 12, 218-220 (2002). [CrossRef]
  22. Y. Rickard and N. Georgieva, "Problem-independent enhancement of PML ABC for the FDTD method," IEEE Trans. Antennas Propag. 51, 3002-3006 (2003). [CrossRef]
  23. S. Cummer, "Perfectly matched layer behavior in negative refractive index materials," IEEE Antennas Wirel. Propag. Lett. 3, 172-175 (2004). [CrossRef]
  24. A. V. Oppenheim and R. W. Schafer, Digital Signal Processing (Prentice-Hall, 1975).
  25. J.-P. Berenger, "Three-dimensional perfectly matched layer for the absorption of electromagnetic waves," J. Comput. Phys. 127, 363-379 (1996). [CrossRef]
  26. O. P. Gandhi, B.-Q. Gao, and J.-Y. Chen, "A frequency-dependent finite-difference time-domain formulation for general dispersive media," IEEE Trans. Microwave Theory Tech. 41, 658-665 (1993). [CrossRef]
  27. W. D. Hurt, "Multiterm Debye dispersion relations for permittivity of muscle," IEEE Trans. Bio. Eng. 32, 60-64 (1985). [CrossRef]
  28. K. E. Oughstun, "Computational methods in ultrafast time-domain optics," Comput. Sci. Eng. 5, 22-32 (2003). [CrossRef]
  29. J. Xi and M. Premaratne, "Analysis of the optical force dependency on beam polarization: Dielectric/metallic spherical shell in a Gaussian beam," J. Opt. Soc. Am. B 26, 973-980 (2009). [CrossRef]
  30. A. Vial, A.-S. Grimault, D. Macías, D. Barchiesi, and M. L. 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(1-7) (2005). [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.

Figures

Fig. 1. Fig. 2. Fig. 3.
 
Fig. 4.
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited