FDTD Modeling of Lorentzian DNG Metamaterials by Auxiliary Differential Equation Method


In this paper, Finite Difference Time Domain (FDTD) is utilized to simulate metamaterials of Double Negative (DNG) origin that refers to those materials having simultaneous negative permittivity and permeability. The problem regarding space formulation is achieved by means of auxiliary differential equation method (ADE), which is easy, reliable and also causal process in nature thus making it proficient. It uses fair approximations to explicate the model. Mur’s boundary condition is used for 1-D problem space and convolution perfectly matched layer boundary is implemented for 2-D problem space. The properties of metamaterial conform their speculations of energy absorption, enhancement and backward propagation property with the aid of graphs engineered by Matlab simulation both in 1-D and 2-D. Also, the interaction of fields on DNG and Double Positive (DPS) layers is contrasted. The results achieved elucidate the validity and effectiveness of the ADE method and the Convolution Perfectly Match Layer (CPML) in designing DNG metamaterials.

Share and Cite:

Goswami, C. , Mukherjee, S. , Karmakar, S. , Pal, M. and Ghatak, R. (2014) FDTD Modeling of Lorentzian DNG Metamaterials by Auxiliary Differential Equation Method. Journal of Electromagnetic Analysis and Applications, 6, 106-114. doi: 10.4236/jemaa.2014.65010.

1. Introduction

Mediums having simultaneously negative permittivity and permeability known as left handed materials (LHM), as conceptualized by Veselago in 1968 [1] , has found numerous applications in microwave circuits and antenna. The LHMs provide some unfamiliar and noteworthy phenomena like negative refraction, reversal of Doppler Effect and Cerenkov radiation etc. This has initiated development of efficient numerical techniques to analyze electromagnetic wave propagation in such left-handed media. The Finite Difference Time Domain (FDTD) method [2] , is an efficient and robust technique which is widely used for modeling electromagnetic wave interaction with various frequency dispersive and non-dispersive materials successfully. The existing frequency dispersive FDTD methods can be broadly classified into three classifications: the Recursive Convolution (RC) method, the Auxiliary Differential Equation (ADE) method and Z-transform method.

In 1990, Luebbers et al. proposed the first frequency dependent FDTD formulation for the modeling of Debye media [3] using the RC method and also used to study wave propagation for Drude materials [4] and Nth order dispersive media [5] . In 2008, D. B. Ge analyzed electromagnetic propagation through dispersive [6] and lossy mediums [7] , implementing FDTD. The bi-isotropic or Chiral media was also investigated using RC-approach [8] in 2004.

Kashiwa and co-workers first used ADE method to formulate Debye media, Lorentz media and Cole-Cole dispersive media [9] -[12] in 1990. Joseph et al. [13] and Gandhi et al. [12] -[14] independently developed an ADE model of Debye media and Nth order dispersive media respectively. A nonlinear phenomenon of dispersive media was investigated by Taflove [15] at the same time.

Z-transforms [16] [17] based on dispersive FDTD formulation was studied first by Sullivan in 1992. Feise et al. [18] proposed Pseudo-Spectral Time Domain (PSTD) technique for the modeling of backward-wave metamaterial and Lee et al. [19] investigated the effective medium approach to modeling split ring resonator LHMs using PLRC method. Ziolkowski also investigated the wave propagation in Drude type [20] and Lorentz type [21] LHMs media in 2001. Chen et al. [22] implemented Z-transform to model Lorentzian DNG metamaterials in FDTD ground in 2004. Then a lot of work has been done in this field. At its very latest, in 2011, K.H. Lee implemented FDTD method on Lorentz-Drude Dispersive Model [23] .

For the truncation of computational domain, the MUR boundary condition and the Perfectly Matched layers (PMLs), invented by Berenger [24] could absorb waves at all frequencies as well as all angles for the linear, anisotropic and other complicated mediums. In 2000, Roden and Gedney invented the CPML (Convolution PML) [25] , which was applicable to complex media including anisotropic media. CPML can also be referred to as Complex frequency shifted PML (CFS-PML) as frequency shifting is accomplished in the argon plane. CPML boundary is also implemented for Lorentzian DNG Metamaterial formulation by Auxiliary Differential Equation Based ADI-FDTD [26] .

In this paper, the formulation of ADE scheme for FDTD formulation in normal DPS as well as DNG based materials is carried out. Section 2 deals with the formulation of ADE method for metamaterials media and followed by boundary condition of MUR and CPML formulation in Section 3; the numerical results of 1-D and 2-D case are discussed in Section 4. The concluding remarks are given in Section 5.

2. ADE Formulation of DNG Metamaterial

2.1. Discrete Maxwell’s Equations and FDTD

The time-dependent Maxwell’s curl equations that are indispensable in FDTD formulation are given in (1) and (2), respectively.



For simplification of the problem, it is assumed that the wave propagates through a vacuum medium. Additionally, the cells of the domain have been assumed to be of the same size. The electromagnetic-field components are arranged on the cells in the same way as that using the conventional FDTD method. Taking the central difference approximations, for both the time and space co-ordinates, we get the FDTD formulations.



where is the space increment and the time step.

The equations shown above allow updating of both the electric and magnetic field in the DPS region; however the formulation changes in the DNG region.



and are such chosen so that at a certain frequency, these parameters become negative, hence modeling the layer DNG.

2.2. The Metamaterial Region

The equation governing the formulation of the metamaterial region are given as in [13]



where and are the free space permittivity and permeability, and are the electric and magnetic plasma frequency, and are the low frequency edge of electric and magnetic forbidden band, and are the electric and magnetic collision frequency.

Now, for calculation of fields in this layer, an inverse Fourier Transform (IFT) to change the s-domain equations into time domain auxiliary differential equations is used. It is further reduced to recursive formulations by the central difference approximation method and hence the fields are updated by time marching algorithm. A new parameter is introduced in the formulation.


where (10)

or (11)

By using the conversion and, and by denoting

as, by recursion in time domain we finally get the updated equations,


Arranging like terms, we get,

or, (13)

Thus, we obtain by this formulation.

Similarly converting (9) into time domain by IFT method we get,


Finally, by the above formula we obtain the electric field in the time domain for the DNG slab region.

3. Boundary Conditions for 1-D and 2-D

To deal with open domain sources, some boundaries are needed as mentioned above to restrict the reflection of the outgoing source waves.

3.1. Absorbing Boundary Conditions (Mur’s ABC)

The ABCs have been used for their easier implementation in the 1-D problem space. At the boundaries, it takes double time steps for a wave front to cross the cell. Since boundaries have been given at both ends, need both upper and lower boundary conditions in order to implement the above situation. The fields are updated accordingly keeping in mind that consistency is maintained throughout the Yee cell. The formulation code used is


3.2. Convolution Perfectly Matched Layer (CPML)

The Mur’s ABC is inefficient and inaccurate for 2-D problem space simulation. In such situations an improve version of Beranger’s PML, the CPML (Convolution PML) is introduced by Roden and Gedney in 2000. CPML is applicable to general media including anisotropic media, while implementing the splitting field technique also. CPML can also be referred to as Complex frequency shifted PML (CFS-PML) as frequency shifting is accomplished in the argon plane. The final electric field components can be expressed as:






These codes are implemented in MATLABTM and simulated to obtain the results for both 1-D and 2-D cases. The code contains various parameters, iterations and matrices that are needed to update the necessitated fields. The simulation was performed in an Intel Core I3 processor, 300 GB Memory, 3 GB RAM using MATLABR2010 version.

4. Numerical Results and Discussions

In this section the 1-D and 2-D numerical results of the above formulation are represented.

4.1. One Dimensional Case

In case of one dimensional, the medium is considered as Figure 1; where the DNG slab is inserted between the two DPS (ε > 0, μ > 0) layer. To formulate the DNG layer, the following data are considered.

The negative unity value arises for permittivity and permeability when. The problem space is considered to be of 300 cells. is chosen as 0.0377, from the Courant condition, where c denotes speed of light.

The DNG slab is interjected from cell no. 120 to cell no. 180, between two DPS slabs each of 120 cells. The rest of the problem space consider as free space. Mur’s Absorbing Boundary condition (ABC), which is very efficient for the case of 1-D problem, is applied at the ends to truncate the problem space.

The Gaussian modulated pulse, whose central frequency equal to ωpeak , parameters spread as 30 and delay as 5, is launched 5 cells from the left boundary of the DPS 1 region.

Comparative study is made between DNG and DPS region where the DNG region of the Figure 1 is replace by DPS (actually it is DPS2 as shown in Figure 2(a)) region with permittivity of 2 which is basically normal dielectric region. The snap shot of electric field distribution at 4995th time step for DPS-DPS-DPS interface and DPS-DNG-DPS interface are shown in Figure 2(a) and Figure 2(b) respectively. A distinguishable comparison from the snap shot is observed for the normal incident of wave. The amplitude of the electric field gradually increases when the wave approach towards the DNG region. It signifies that the DNG medium has a capability to increase the energy of the wave compared to normal DPS region.

Figure 1. Geometry of 1-D DPS-DNG-DPS problem space where wave propagation along x direction.


Figure 2. Wave electric field distribution at time step t = 4995th (a) the DPS-DPS-DPS interface (b) the DPS-DNG- DPS interface.

The time histories of the electric field recorded at two specific points within the DNG medium, one at the center of the DNG medium (x = 150∆x) and another at 20 cells ahead (x = 130∆x) from the center. The electric field profile and are recorded for early time and delay time as shown in Figure 3(a) and Figure 3(b) respectively. It is observed that the early time response at point 1 is occurs before the point 2. This shows that causality in the direction of the wave propagation is preserved in the DNG medium, which is in agreement with the results presented in [20] . Another interesting property observed in Figure 3 is that the response at the second point leads the one at the first point in the delay time

4.2. Two Dimensional Case


Figure 3. Time histories of FDTD electric field at the x = 150∆x center of the DNG slab and at x = 130∆x, 20 cells ahead from the center (a) for early-time results; (b) for delay-time results.

Figure 4. Geometry of DNG embedded 2-D problem space implementing CPML.

The interaction of TM wave with DNG slab is display in Figure 5 at different time instants. The wave is launch from DPS region, 3 cells before the CPML boundary. The wave is incident normally to the DNG slab and after entering the DNG region, the backward wave propagation is observed, which is one of the characteristics of left-handed metamaterial. Again forward propagation of the wave is observed when it reaches at the second DPS region.

5. Conclusion

The FDTD formulations incorporated with auxiliary differential equation method is presented for modeling the Lorentzian DNG metamaterial. The effectiveness of ADE method is established for the case of complex media. Mainly two properties of the metamaterial are discussed in 1-D and 2-D FDTD formution. The DNG metamaterial has capability to increase the amplitude of field component, compared to normal dielectric media, when the wave is passed through it which is observed from 1-D simulations. The causality is also maintained throughout the metamaterial layer and the energy gets absorbed in DNG medium from the surrounding media and hence behaves as a secondary energy source. The 2-D results demonstrate the backward propagation of wave in DNG media. An improvement over the conventional PML, CPML has been used here as efficient boundary conditions. As, the PML is inefficient in absorbing evanescent waves, the source must be placed faraway from the

(a) (b)(a) (b)

Figure 5. Electric field distribution at different time step of 2-D DPS-DNG DPS interface (a) at t = 75∆t, (b) at t = 250∆t, (c) at t = 359∆t, (d) at t = 855∆t.

boundary. This increases the number of cells and hence computational requirements and memory. Our results demonstrate that CPML performs very well even with the source being introduced close to the boundaries. Also PML recoreds late time reflections due to its weak causal nature. Specially, the 2-D implemention for a metamaterial is very difficult for PML as PML takes into account a continous function all throughout the region. CPML being strictly causal has no problem regarding evanescent waves and hence, source is placed very close to the boundaries, thus saving memory [8] .

Conflicts of Interest

The authors declare no conflicts of interest.


[1] Veselago, V.G. (1968) The Electrodynamics of Substance with Simultaneously Negative Values of ε and μ. Soviet Physics Uspekhi, 10, 509. http://dx.doi.org/10.1070/PU1968v010n04ABEH003699
[2] Taflove, A. (1995) Computational Electrodynamics: The Finite-Difference Time-Domain Method. Artech House, Norwood.
[3] Luebbers, R.J., Hunsberger, F., Kunz, K.S., Standler, R.B. and Schneider, M. (1990) A Frequency Dependent Finite-Difference Time-Domain Formulation for Dispersive Materials. IEEE Transactions on Electromagnetic Compatibility, 32, 222-227. http://dx.doi.org/10.1109/15.57116
[4] Luebbers, R.J., Hunsberger, F. and Kunz, K.S. (1991) A Frequency-Dependent Finite-Difference Time-Domain Formulation for Transient Propagation in Plasma. IEEE Transactions on Antennas and Propagation, 39, 29-34.
[5] Luebbers, R.J. and Hunsberger, F. (1992) FDTD for Nth-Order Dispersive Media. IEEE Transactions on Antennas and Propagation, 40, 1297-1301. http://dx.doi.org/10.1109/8.202707
[6] Zhang, Y.-Q. and Ge, D.-B. (2009) A Unified FDTD Approach for Electromagnetic Analysis of Dispersive Objects. Progress in Electromagnetics Research, PIER96, 155-172.
[7] Jiang, Y.-N., Ge, D.-B. and Ding, S.-J. (2008) Analysis of TF-SF Boundary for 2D-FDTD with Plane P-Wave Propagation in Layered Dispersive and Lossy Media. Progress in Electromagnetics Research, PIER83, 157-172.
[8] Akyurtlu, A. and Werner, D.H. (1990) A Novel Dispersive FDTD Formulation for Modeling Transient Propagation in Chiral Metamaterials. IEEE Transactions on Antennas and Propagation, 52, 2267-2276.
[9] Kashiwa, T., Yoshida, N. and Fukai, I. (1990) A treatment by the Finite-Difference Time-Domain Method of the Dispersive Characteristics Associated with Orientation Polarization. IEEE Transactions on Antennas and Propagation, E73, 1326-1328.
[10] Kashiwa, T. and Fukai, I. (1990) A Treatment by the FD-TD Method of the Dispersive Characteristics Associated with Electronic Polarization. Microwave and Optical Technology Letters, 3, 203-205.
[11] Kashiwa, T., Ohtomo, Y. and Fukai, I. (1990) A Finite-Difference Time-Domain Formulation for Transient Propagation in Dispersive Media Associated with Cole-Cole’s Circular ARC Law. Microwave and Optical Technology Letters, 3, 416-419. http://dx.doi.org/10.1002/mop.4650031204
[12] Joseph, R.M., Hagness, S.C. and Taflove, A. (1991) Direct Time Integration of Maxwell’s Equations in Linear Dispersive Media with Absorption for Scattering and Propagation of Femtosecond Electrogmagnetic Pulses. Optics Letters, 16, 1412-1414. http://dx.doi.org/10.1364/OL.16.001412
[13] Gandhi, P., Gao, B.Q. and Chen, J.Y. (1992) A Frequency-Dependent Finite-Difference Time Domain Formulation for Induced Current Calculations in Human Beings. Bioelectromagentics, 13, 543-556.
[14] Gandhi, P., Gao, B.Q. and Chen, J.Y. (1993) A Frequency-Dependent Finite-Difference Time Domain Formulation for General Dispersive Media. IEEE Transactions on Microwave Theory and Techniques, 41, 658-665.
[15] Goorjian, P.M. and Taflove, A. (1992) Direct Time Integration of Maxwell’s Equations in Nonlinear Dispersive Media for Propagation and Scattering of Femto Second Electromagnetic Solutions. Optics Letters, 17, 180-182.
[16] Sullivan, D.M. (1992) Frequency-Dependent FDTD Methods Using Z Transforms. IEEE Transactions on Antennas and Propagations, 40, 1223-1230. http://dx.doi.org/10.1109/8.182455
[17] Demir, V., Elsherbeni, A.Z. and Arvas, E. (2005) FDTD Formulation for Dispersive Chiral Media Using the Z Transform Method. IEEE Transactions on Antennas and Propagations, 53, 3374-3384.
[18] Feise, M.W., Schneider, J.B. and Bevelacqua, P.J. (2004) Finite-Difference and Pseudospectraltime-Domain Methods Applied to Backward-Wave Metamaterials. IEEE Transactions on Antennas and Propagations, 52, 2955-2962.
[19] Lee, J.Y., Lee, J.H., Kim, H.S., Kang, N.W. and Jung, H.K. (2005) Effective Medium Approach of Left Handed Material Using a Dispersive FDTD Method. IEEE Transactions on Magnetics, 41, 1484-1487.
[20] Ziolkowski, R.W. and Heyman, E. (2001) Wave Propagation in Media Having Negative Permittivity and Permeability. Physical Review E, 64, Article No. 056625. http://dx.doi.org/10.1103/PhysRevE.64.056625
[21] Ziolkowski, R.W. (2001) Superluminal Transmission of Information through an Electromagnetic Metamaterial. Physical Review E, 63, Article No. 046604.
[22] Suwailam, M.M.B. and Chen, Z.Z. (2004) FDTD Modeling of Lorentzian DNG Meta-Materials with the Z-Transform. In: 3rd International Conference on Computational Electromagnetics and Its Applications Proceedings, IEEE, New York, 40-43.
[23] Lee, K.H., Ahmed, I., Goh, R.S.M., Khoo, E.H., Li, E.P. and Hung, T.G.G. (2011) Implementation of the FDTD Method Based on Lorentz-Drude Dispersive Model on GPU for Plasmonic Applications. Progress in Electromagnetics Research, 116, 441-456.
[24] Berenger, J.P. (1994) A Perfectly Matched Layer for the Absorption of Electromagnetic Waves. Journal of Computational Physics, 114, 185-200. http://dx.doi.org/10.1006/jcph.1994.1159
[25] Roden, A. and Gedney, S.D. (2000) Convolutional PML (CPML): An Efficient FDTD Implementation of the CFSPML for Arbitrary Media. Microwave and Optical Technical Letter, 27, 334-339.
[26] Mukherjee, S., Karmakar, S., Goswami, C. and Ghatak, R. (2012) Electromagnetic Wave Propagation Modeling in Lorentzian DNG Metamaterial by Auxiliary Differential Equation Based ADI-FDTD. In: Proceedings of 2012 1st International Conference on Emerging Technology Trends in Electronics, Communication and Networking, IEEE, New York, 1-4.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.