## Abstract

We present an auxiliary differential equation Finite-difference Time-domain (ADE-FDTD) approach to numerically model the wave propagation within a gain or absorbing medium such as quantum well structures. Start from traditional quantum electronics theory, the macroscopic susceptibility of the semiconductor is derived and expressed by a multiple-Lorentz-like model based on Prony’s method. With the auxiliary differential equation method each Lorentz-like model can be simulated in the time domain and the induced polarization is then determined by summing all the models. By incorporating the induced polarization into the time-domain Maxwell’s equations, electromagnetic wave propagation in the quantum well medium can be accurately modeled using the FDTD method.

©2006 Optical Society of America

## 1. Introduction

During the last few decades, semiconductor devices have triggered a significant level of interest, due to the wide variety of applications in which they are used, such as systems associated with communications, remote sensing, and optical computing [1–3]. Most current simulation models for such devices lack the ability to accurately simulate the physics within a gain or absorbing region since they are based on the known macroscopic and phenomenological models. Recently, Amit and Robert^{4} and Ziolkowski^{5} have accurately simulated the electromagnetic interactions in a dispersive and nonlinear medium by incorporating macroscopic polarization into Maxwell’s equations, where the quantum mechanical effects and time-domain dynamics of the atomic populations are governed by rate equations. In fact, their two-level or four-level atom models can be considered as Lorentz models. However, for most semiconductor materials, the spectral broadening of transition of the photons has to be taken into account. In this case the single Lorentz model is no longer valid and an alternative model must be used.

In this paper we present such a method and begin with the semi-classical density matrix theory [6, 7] to derive the quantum mechanical descriptions of the semiconductor medium. Such physical models incorporate all propagation effects so that an associated macroscopic susceptibility can be derived. We then couple Maxwell’s equations with the induced polarization to numerically solve the quantum well medium. Usually, the susceptibility of quantum well is very complex due to the complicated conduction and valence band structures, which bring about the invalidity of the single Lorentz model in such materials. In addition, when the injection current increases, the response of the susceptibility and change are no longer linear. So, the key point is to accurately describe the susceptibility. In our model, we use the Fourier transformation to derive the time response of the susceptibility. Then the Prony’s method [8] is applied to approximate it with a superposition of exponential series so that the susceptibility of a quantum well can be expressed as a superposition of multiple Lorentz-like models. Using the inverse Fourier transform one can derive the time-domain response. In this way, each term of multiple Lorentz-like models can be solved using FDTD technique, as we show.

Accordingly, the remainder of the paper is organized as follows: In Sec. II, based on the semi-classical density-matrix theory, the microscopic susceptibility in semiconductor structure is obtained. Then, multiple Lorentz-like models are derived using the Prony’s method. In Sec. III, the FDTD method is applied to numerically solve the quantum well structure, and the numerical results are provided in Sec. IV. Finally, a conclusion and summary are presented in Sec. V.

## 2. Macroscopic susceptibility

Based on the semi-classical density matrix theory, we consider a two-level system of semiconductor medium with a rigorous **k** selection rule applied to the recombining electron-hole pairs, the band structure of which is shown in Fig. 1. In this case, the density-matrix equations can be written in terms of the distribution (occupation) functions for electrons *f*_{e}
and holes *f*_{h}
, in the presence of an optical field *E*(*t*), as

where “^{*}” stands for the complex conjugate, and

and

are the quasi-Fermi distributions that the electrons and the holes tend to relax to in the presence of the optical field perturbation, *ħ* is Plank’s constant /2*π* , *τ*_{e}
and *τ*_{h}
are the intraband relaxation time of the electrons in conduction band and holes in valence band, respectively. *E*_{e}
(*k*) and *E*_{h}
(*k*) are the energies of the electron and hole with wave vector *k*, respectively. *F*_{e}
and *F*_{h}
are quasi-Fermi energy levels for electrons and holes, respectively, and *ρ*_{eh}
(*k*) is off-diagonal (electron-hole) density matrix element, *E*_{k}
is possible transition energy of a electron-hole pair with wave vector *k*, *T*
_{2} is the intraband dephasing time and *γ* is the dipole moment matrix element. ${f}_{e}^{0}$(*k*) and ${f}_{h}^{0}$(*k*) are implicitly dependent on the optical and injection current in the semiconductor region. The steady-state solution of Eqs. (1) and (2) can be written as

where *L*(*ħω*-*E*_{k}
) = *E*
^{2}
_{T2} /[*E*
^{2}
_{T2} + (*ħω* - *E*_{k}
)^{2}] is the Lorentz (broadening) function, *E*
_{T2} = *ħ*/*T*
_{2} is the space average photon density inside the active region, *I*_{s}
= ħ^{2}
*ε*
_{0}
*n*
^{2}/*γ*
^{2}
*E*(*τ*_{e}
+ *τ*_{h}
)*T*
_{2} is the saturation photon density, and *n* is the model refractive index. *f*_{e}
(*k*) and *f*_{h}
(*k*) are the nonlinear function of the presence field or average photon density. Under the small signal approximation, in which the photon density is very small compared with saturation density, the second terms in Eq. (6) and Eq. (7) can be omitted, thus *f*_{e}
(*k*) and *f*_{h}
(*k*) can be simplified to ${f}_{e}^{0}$(*k*) and ${f}_{h}^{0}$(*k*), respectively. In this case, by performing the Fourier transform, the solution of Eq. (3) is

where “⊗” stands for convolution. As is well known, the induced polarization can be written quantum mechanically as

or classically as

where *χ* is the macroscopic susceptibility. From Eq. (9) and Eq. (10), we obtain the expression in frequency domain as

Combine Eq. (11) with Eq. (6)–Eq. (8), the susceptibility can be derived as

In this case, *V* stands for the volume of the active medium region. Please be advised that this is valid for small signal intensities and slowly varying pulses.

Next, we consider a quasi-two-dimensional (2D) system for the movement of the charge carriers in the quantum well (QW) semiconductor. [9] In the horizontal plane (*x, y* plane) of the QW, the movement remains unrestricted; in the vertical dimension *z*, it is confined. The density of states in such a system is obtained from the energy levels *E*_{n}
(*n* = 1,2,⋯) of the quantized dimension *z*, each forming a subband *n* with respect to the free *x, y* movement. The 2D density of states *D*(*E*'} , i.e., the available states per unit area and unit energy, can be written as the sum over the step-like density of states functions of each individual subband. [10]

where *H* is the Heaviside unit step function, *m*_{n}
is the effective mass in the *n*
^{th} 2D sub-band
with respect to the *x, y* movement, *d*_{z}
is the thickness of quantum well structure. The energy *E*' is measured from the band edge of the bulk semiconductor into the band. The quasi-Fermi levels of the electrons and holes can represent the injection level, and, the quasi-Fermi levels for electrons and holes are related through the charge neutrality condition. In the valence band, there are usually two branch holes, namely, the heavy (h) and light (l) holes. As a result of the conservation of in-plane momentum *k*_{x,y}
and confinement quantum number *n* during a dipole transition of injected electrons from conduction to the valence bands, the 2D density of states functions of the two sub-bands with same n can be combined to form a reduced 2D density of states describing the dipoles of energy *E*' = *ħω* between the electron and the respective hole; i.e.,

where the effective mass *m*_{n}
in Eq. (20) is replaced by reduced effective mass,

and

is the gap between the subband pair *e*_{n}
~ *ih*_{n}
. The linear susceptibility *χ* can be rewritten as

The square of transition matrix elements in the QW structures generally can be written as [11]

where *δi* is defined as the polarization modification factor, and ${M}_{b}^{2}$ is the square of the transition matrix element given by :

in which *ξ* is the Matrix element prefactor [12], for GaAs *ξ* is 2.66. The polarization modification factor in the QW provides an enhancement of the oscillator strength for TE polarization at photon energies near the gain peak as opposed to TM polarization, where the oscillator strength is diminished. Thus, for a QW structure, stability of lasing in the TE mode is improved further, and TM polarization need not be considered. The TE polarization modification factors for the e-hh and e-lh transitions, respectively, are:

With the modified matrix element, the susceptibility eventually becomes

The energy levels in the quantum well structures in the *z* direction can be obtained by solving the Schrodinger equation for a one-dimensional potential well, the following eigenvalue equations are obtained [11],

where *V*
_{0} is the depth of the potential well and *m* is the effective mass. The eigenvalue equation can be numerically solved to yield the energy levels *E*_{n}
in a potential well.

## 3. Auxiliary differential equation Finite-Difference Time-Domain (ADE-FDTD) method

In the last section, we obtain the macroscopic susceptibility of semiconductor material based on the semi-classical density matrix theory, thus, the induced polarization can be incorporated into the Maxwell’s equations. In general, the material is dispersive and nonlinear due to the transitions and combinations of electrons and holes. Under the small signal approximation the material becomes linear, as a result we only need consider it as a dispersive material. Once the induced polarization of such a material is determined, we can then apply numerical techniques to simulate the wave propagation in the corresponding gain or absorbing medium. To further simplify the algorithm, we instead use induced polarization, where we derive an equivalent conductivity for the quantum well medium to incorporate into the Maxwell’s equations.

As is well known, the FDTD method has become one of the more widely used numerical techniques for solving electromagnetic boundary value problems due to its reduced memory requirements and its simplicity of implementation. In short, the FDTD method directly discretize Maxwell’s equations in both time and space by using central difference technique, and iteratively updates the electrical and magnetic field components at all points in the computational space using the leapfrog fashion. The material parameters of the medium are also incorporated in the computation. To this end, the method is able to easily deal with the interactions between electromagnetic fields and objects having complex shape and/or inhomogeneous media.

Over the past several years, a number of different methods have been proposed to account for material dispersion with absorption and gain using the FDTD method. These include recursive convolution (RC) methods [13, 14] that implement a discretized convolution of the dispersion relation, ADE methods that discretize a differential equation obtained from the relevant constitutive relations, and *Z*-transform methods [15] based on digital processing techniques. Although the RC method can be easily implemented, it needs the entire field history to implement the convolution, and the convolution involves a very demanding computer burden. As a result, we choose an alternative method, namely the ADE method to solve the problem.

In general, Maxwell’s equations within the quantum well medium can be written in the form

where **P** is the induced polarization. For simplicity, we considered here only a one-dimensional problem with field components *E*_{x}
and *H*_{y}
propagating along *z* direction through a nonmagnetic, isotropic medium. In this case, the above equations reduced to:

where the polarization *P*(*t*) is defined in Eq. (10). Also, the susceptibility *χ*(*t*)in the time domain can be obtained by performing the Fourier transform on Eq. (17), which is given by:

where *ω*_{T}
= 2*π*/*T*
_{2} and *U*(*t*) is the unit step function. From Eq. (25), since *χ*(*t*) is implicit function of *t*, the derivative of *P*(*t*) with respect to *t* can be performed analytically. To this end, we introduce the polarization current density *J*(*t*) to replace the *P*(*t*),

The corresponding frequency field form is *J*(*ω*) = σ(*ω*)*E*(*ω*). The equivalent conductivity σ(*t*) can be derived from $\sigma \left(t\right)=\frac{\partial \chi \left(t\right)}{\partial t}$ as,

Obviously, the response of conductivity σ(*t*) is very complex due to the presence of the quantum well. One reason for this is that, at the higher excitation of the quantum well medium, the transitions involving the second subband can be very important, so that a spectrum with two peaks may appear. Another reason is that the semiconductor medium is constructed by collecting two-level atoms with a range of transition frequencies, which is similar to an inhomogeneous broadened two-level system. As a result, it involves multiple oscillators.

Once we have known the time domain response of σ(*t*), the ADE method is applied. To
do so, we reformulate the frequency dependent conductivity function. It is well known that if the time domain σ(*t*) can be approximated by a sum of complex exponential functions, one can use the highly efficient auxiliary differential equation. Thus, by using Prony’s method it is possible to set up the complex exponentials:

in which *α*_{i}
is damping factor of the *i*-th mode, *ω*_{i}
is the resonant frequency of the *i*-th mode, and *P* is the order of the model. Usually there are polarization pole pairs, therefore, in frequency domain,

where *A*_{i}
+ *jB*_{i}
= *C*_{i}
. Using Prony’s method, we record *N* equally spaced time samples of σ(*N*) in the time domain and set up the following N-p equations:

$${\sigma}_{p+1}+{E}_{1}{\sigma}_{p}+{E}_{2}{\sigma}_{p-1}+\cdots +{E}_{p}{\sigma}_{2}=0$$

$$\vdots \phantom{\rule{2.2em}{0ex}}\vdots \phantom{\rule{2.2em}{0ex}}\vdots \phantom{\rule{2.2em}{0ex}}\ddots \phantom{\rule{2.2em}{0ex}}\vdots \phantom{\rule{2.2em}{0ex}}\vdots $$

$${\sigma}_{N}+{E}_{1}{\sigma}_{N-1}+{E}_{2}{\sigma}_{N-2}+\cdots +{E}_{p}{\sigma}_{N-p}=0$$

Solving this set of equations for *E*_{l}
permits the *D*_{l}
to be found as the roots of the polynomial

The next step is to find the *C*_{l}
in Eq. (28). This is accomplished by writing out Eq. (28) for each of the time samples, thereby obtaining a set of *N* equations:

$${D}_{1}{C}_{1}\phantom{\rule{.2em}{0ex}}+{D}_{2}{C}_{2}\phantom{\rule{.2em}{0ex}}+\cdots +{D}_{p}{C}_{p}={\sigma}_{2}\phantom{\rule{5.2em}{0ex}}$$

$${\left({D}_{1}\right)}^{2}{C}_{1}+{\left({D}_{2}\right)}^{2}{C}_{2}\phantom{\rule{.2em}{0ex}}+\cdots +{\left({D}_{p}\right)}^{2}{C}_{p}+{\sigma}_{3}\phantom{\rule{3.2em}{0ex}}$$

$$\phantom{\rule{.2em}{0ex}}\vdots \phantom{\rule{3.2em}{0ex}}\vdots \phantom{\rule{3.2em}{0ex}}\ddots \phantom{\rule{2.2em}{0ex}}\vdots \phantom{\rule{3.2em}{0ex}}\vdots $$

$${\left({D}_{1}\right)}^{N-1}{C}_{1}+{\left({D}_{2}\right)}^{N-2}{C}_{2}+\cdots +{\left({D}_{p}\right)}^{N-1}{C}_{p}={\sigma}_{N}$$

Once *C*_{l}
and *D*_{l}
are solved, the *α*_{i}
, *ω*_{i}
, *A*_{i}
and *B*_{i}
can be found. Then the polarization electric current density can be expressed as the superposition of multiple polarization current components

By performing an inverse Fourier transform of each term in the above equation, we have

The above equations are required ADEs for the equivalent induced current *J* (*t*). This can be easily and accurately implemented using the FDTD method. The electromagnetic field components and induced electric currents are shown in the grid as indicated in Fig. 2.

Using the central difference expression, we can discretize Eq. (24) and Eq. (35), in which the electric field is centered at time-step *n*+1/2, and magnetic field and polarization current are centered at time-step *n*. They are given by

$$-\frac{\left({A}_{i}{\alpha}_{i}+{B}_{i}{\omega}_{i}\right)\Delta {t}^{2}-2{A}_{i}\Delta t}{1-{\alpha}_{i}\phantom{\rule{.2em}{0ex}}\Delta t}{E}_{x}^{\frac{n+1}{2}}\left(I\right)-\frac{\left({A}_{i}{\alpha}_{i}+{B}_{i}{\omega}_{i}\right)\Delta {t}^{2}+2{A}_{i}\Delta t}{1-{\alpha}_{i}\phantom{\rule{.2em}{0ex}}\Delta t}{E}_{x}^{\frac{n-1}{2}}\left(I\right)$$

Thus, ADE-FDTD algorithm can be used to model a dispersive medium with *P* pairs Lorentz pole. The reformulated ADE method eliminates the requirement to solve for the field history, and, it needs a smaller number of unknowns to store than the RC method. In the following section, we present a few of examples to illustrate the utility of the proposed algorithm.

## 4. Numerical results

We first present, using Prony’s method to numerically solve the parameters in Eq. (29), a comparison simulation result of the quantum well structure based on our ADE-FDTD method with its analytical result. Following this, a more general application is given.

#### 4.1 Numerically solving macroscopic susceptibility with Prony’s method

Our numerical calculations are performed for a representative separate confinement quantum well structure in the GaAs/AlGaAs system. The parameters for the material used in the calculations are listed in Table 1 as follows:

Figure 3 shows a comparison between the results using Prony’s method, with finite number of pairs of Lorentz-like models, and the analytical solution of susceptibility *χ*(*ω*). In our simulation, we set *P* in Eq. (29) as 20, which means that superposition of 10 pairs of Lorentz-like models are used to approximate the analytical *χ*(*ω*). The injection current density *N* is equal to 6×10^{18}/ cm^{-3}. As shown in Fig. 3, Prony’s method data shows extremely good agreement with the analytical result based on quantum electronics theory given in Eq. (21). Therefore, the obtained parameters in the Lorentz-like model can be used for the FDTD simulations.

#### 4.2 ADE-FDTD simulation

By using Prony’s method, the dispersive property of a quantum well structure, such as equivalent conductivity, can be approximated by the superposition of multiple Lorentz-like models. The modified ADE-FDTD is then applied to simulate the wave propagation in a homogeneous quantum well material to examine the presented scheme. To do this, we consider a sinusoidal modulated Gaussian pulse with the center frequency *f*_{c}
of 4.5 × 10^{14} Hz and a homogeneously broadened bandwidth Δ*f*_{c}
of 2.5×10^{14} Hz propagating through the quantum well structure. The FDTD grid size is $\Delta =\frac{c}{200\left({f}_{c}+\Delta {f}_{c}\right)n}=0.6\mathit{nm}$ in which the refractive index *n* is equal to 3.6. The time step Δ*t* is chosen as 1.8×10^{-18}s so that the time step is slightly smaller than the value dictated by the Courant condition. The electromagnetic field was recorded at two different locations that were 100 grid points apart. Then, the recorded temporal fields were Fourier transformed to calculate the frequency spectrum of the medium. Figures 4 shows a comparison between ADE-FDTD and analytical result with N = 1×10^{18}/ cm^{-3}, 6×10^{18}/ cm^{-3} and 12×10^{18}/ cm^{-3}, respectively. As shown the FDTD results show good agreement with the analytical results. Also, the figures show that as the injection current density increase, the gain character of the absorption factor becomes increasingly clearer. When *N* = 1×10^{18}/ cm^{-3}, there is no gain effect, while with *N* increasing to 6×10^{18}/ cm^{-3} and 12×10^{18}/ cm^{-3}, the peak of the absorption factor also increases to 1.02 and 1.04, respectively.

#### 4.3. Application

In this section, we consider the electromagnetic wave propagation of a pulse through a multiquantum well (MQW) structure, as shown in Fig. 5(a). The structure consists of three parts. Part (1) and (3) are homogeneous 35% AlGaAs material without QWs while part (2) is formed by symmetrically imbedding GaAs QWs in the barrier material of 35% AlGaAs. The spacing between the adjacent quantum well is set as 30nm which is equal to the thickness of the quantum well so that the quantum wells are uniformly distributed in part (2). We use 40 layers of quantum wells, which mean that the length of part (2) is 2.4*μm*. Parts (1) and (3) are also set to be 2.4*μm* so the total length of our structure is 7.2*μm*.

To model this, we launch a sinusoidal modulated Gaussian pulse with central frequency located at 450THz and bandwidth equal to 250THz, which is the same frequency range as what we used to calculate the parameters in ADE-FDTD algorithm. The injection current density of the quantum well is set to be *N*= 6×10^{18}/ cm^{-3}. To minimize the reflection on the boundary, the perfectly matched layer absorbing boundary condition is used.[16] In the simulation, the snapshot of the electric field is recorded every 0.015ps. Figure 5(b) shows the evolution of the pulse as it propagates along the medium. The region between two vertical dashed lines indicates the MQW structure. Since part (1) is homogeneous 35% AlGaAs material, the pulse shape remains unchanged in the first two snapshots. Once the pulse enters part (2), which is inhomogeneous MQW medium, the pulse shape begins to become broadened as depicted in Fig. 5(b). Such broadening phenomenon happens due to the finite lifetimes of a given energy state[11]. If we assume the state decays exponentially with time, then the Fourier energy spectrum necessary to construct the time-dependent state has a Lorentzian lineshape; hence the energy of each state is no longer sharp but has energy spread over a range of Δ*E*. The Lorentzian lineshape function or broadening function has been introduced in Eq. (6) and Eq. (7). Please be advised that the tiny wavelet on the long tail of the pulse (especially clear in snapshot number 6) is due to the multiple reflections on the boundary of single quantum well layer and adjacent barrier 35% AlGaAs medium. Recall that we have a sandwich-like structure in part (2); each sandwich layer will form a cavity since the material on each side of the small cavity is different from that inside the cavity region. Although the refractive indices are all equal to 3.6, which is the same as the refractive index of GaAs, the permittivity in quantum well layer has an imaginary part while in the adjacent homogeneous 35% AlGaAs layer has not. This is the reason of the formation of those tiny wavelets. This pulse tail may be minimized if we consider a narrower frequency range, for instance, only the frequency range that offers gain. After part (2), the shape remains the same since part (3) is also homogeneous similar with part (1).

By collecting the temporal data of two detectors 6*μm* apart located in part (1) and part (3), respectively, we can calculate the frequency response of this structure. As shown by Fig. 5(c), the absorption factor is quite similar with Fig. 5 while the phase shift is different. This is easy to understand since we have a non-homogenous material and the sandwich-like structure may change the phase. For the high frequency range of the absorption factor, say from 600THz to 700THz, the curve is not as flat as that in Fig. 5. This is due to the cavity effect as explained above. Since the cavity dimension is quite small, the high frequency electromagnetic waves tend to build up a standing wave, or cavity mode, which is the reason that the oscillation of the absorption is more obvious near 700THz.

## 5. Summary

In summary, we combine quantum electronics theory with the ADE-FDTD method to numerically simulate semiconductor devices. The two-level atomic model includes the spectral broadening of transition of the photon. Prony’s method is used to characterize the macroscopic expression of the quantum well structure by a superposition of decaying exponentials. The resultant multiple Lorentz model forms auxiliary differential equations for equivalent electrical currents. The FDTD method is then applied to accurately model the quantum well structures. The results show very good agreement with current theoretical models for the small signal frequency response. Although the method is an approximation under the small signal condition, by properly rearranging Eq. (6) and Eq. (7), it can also be generalized to overcome such a limit. We are currently working on this issue.

## References and links

**1. **H. P. Sardesai and A. P. Weiner, “Nonlinear fiber-optic receiver for ultra short pulse code division multiple access communications,” IEEE Electron. Lett. **33**, 610–611 (1997). [CrossRef]

**2. **C. A. Kapetanakos, B. Hafizi, H. M. Milehberg, P. Sprangle, R. F. Hubbard, and A. Ting, “Generation of high-average-power ultrabroad-band infrared pulses,” IEEE J. Quantum Electron. **35**, 565–576 (1999). [CrossRef]

**3. **L. E. M. Brackenbury, “Multiplexer as a universal computing element for electro-optic logic systems,” IEEE Proc. J. Optoelectron. **137**, 305–310 (1990). [CrossRef]

**4. **A. S. Nagra and R. A. York, “FDTD analysis of wave propagation in nonlinear absorbing and gain media,” IEEE Tran. Antennas Propag. **46**, 334–340 (1998). [CrossRef]

**5. **R. W. Ziolkowski, J. M. Arnold, and D. M. Gogny, “Ultrafast pulse interactions with two-level atoms,” Phys. Rev. A **52**, 3082–3094 (1995). [CrossRef] [PubMed]

**6. **M. Sargent, M. O. Scully, and W. E. lamb, *Laser Physics* (Addison-Wesley, Reading, MA, 1974).

**7. **A. Yariv, *Quantum Electronics*, 3rd ed. (Wiley, New York, 1989).

**8. **F. B. Hildebrand, *Introduction to Numerical Analysis* (Dover, New York, 1974).

**9. **D. Kasemset, C. S. Hong, N. B. Patel, and P. D. Dapkus, “Graded barrier single quantum well lasers-theory and experiment,” IEEE J. Quantum Electron. **19**, 1025–1030 (1983). [CrossRef]

**10. **R. Dingle and I. Festkorperprobleme, *Advances in solid state physics* (Braunschweig, Pergamon-Vieweg).

**11. **P. S. Zory, “Quantum Well Lasers,” (Academic Press INC, CA, 1993).

**12. **N. K. Dutta, “Calculated threshold current of GaAs quantum well lasers.,” J. Appl. Phy. **53**, 7211–7214 (1983). [CrossRef]

**13. **R. J. Luebbers and F. Hunsberger, “FDTD for Nth-order dispersive media,” IEEE Trans. Antennas Propag. **40**, 1279–1301 (1992). [CrossRef]

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

**15. **D. M. Sullivan, “Frequency-dependent FDTD methods using Z transforms,” IEEE Trans. Antennas Propag. **40**, 1223–1230 (1992). [CrossRef]

**16. **J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**, 185–200 (1994). [CrossRef]