Symplectic FDTD method study left-handed material electromagnetic characteristics

Optik 124 (2013) 1716–1720

Contents lists available at SciVerse ScienceDirect

Optik journal homepage: www.elsevier.de/ijleo

Symplectic FDTD method study left-handed material electromagnetic characteristics Hong Wei Yang ∗ , Haiyan Song Department of Physics, Nanjing Agricultural University, Nanjing 210095, People’s Republic of China

a r t i c l e

i n f o

Article history: Received 9 December 2011 Accepted 12 May 2012

Keywords: Symplectic propagation technique Symplectic finite-difference time-domain (SFDTD) Left-handed materials Drude

a b s t r a c t A symplectic finite-difference time-domain (SFDTD) method is used. The method is based on expressing the dynamical system in Hamiltonian form, it can maintain the energy conservation (i.e. Hamiton volume conservation), and the corresponding mathematical formulation for the evolution of the system state is symplectic transformation. This method has better stability and accuracy than the general FDTD method. This is the first use of the method to study the electromagnetic properties of left-handed materials. Using the Drude dispersion medium model to describe the dispersion characteristics of left-handed materials, and left-handed materials with negative refractive index, the collection effect and phase compensation effect are simulated. © 2012 Elsevier GmbH. All rights reserved.

1. Introduction The electrical permittivity ε and the magnetic permeability  are the fundamental characteristic quantities that determine the propagation of electromagnetic waves in matter. Left-handed materials with unusual electromagnetic properties, such as double negative (DNG) feature, the left hand features, backward wave property, negative refraction and the reversed Doppler effect and the reversed Cerenkov radiation effects [1,2]. Symplectic finite-difference time-domain (SFDTD) method is based on dynamic system expressed in Hamilton, it can maintain the energy conservation (i.e. Hamilton volume conservation), and the corresponding mathematical formulation for the evolution of the system state is symplectic transformation. SFDTD method considers Maxwell’s equations as an infinite-dimensional Hamilton system. The Hamilton method should be produced within the symplectic geometry frame. Evolved with time, the system is always a symplectic transformation evolution, so it keeps the Maxwell’s equations symplectic structure. Therefore, this method also has good stability and high accuracy, and compared with the FDTD [3–6] methods for validation. In this paper, Drude model [7–9] as a model of electromagnetic parameters of left-handed materials, the first application of SFDTD algorithms on propagation characteristics of the electromagnetic wave through left-handed materials is analyzed theoretically, and the propagation process is simulated, so as to prove the negative refraction effect, a clustering

effect and the effect of phase compensation of the left-handed materials. 2. SFDTD method theory In case of lossless media, canonical form of Maxwell’s equations is as follows: ∂H ∂t

=−

∂E ∂t

=

0030-4026/$ – see front matter © 2012 Elsevier GmbH. All rights reserved. http://dx.doi.org/10.1016/j.ijleo.2012.06.009

(1)

1 ∇ ×H ε

Maxwell’s equations in electromagnetic field can be expressed with the Hamilton function. Through the variational method (1) can be rewritten as ∂ ∂t

  E



U=

 

= (U + V )

H

E

(2a)

H

{0}3×3

−−1 {R}3×3

{0}3×3

{0}3×3

⎛ 0

⎜ ⎜ ⎜ ∂ R=⎜ ⎜ ∂z ⎝ −

∗ Corresponding author. E-mail address: phd [email protected] (H.W. Yang).

1 ∇ ×E 

∂ ∂y

−

∂ ∂z

∂ ∂y

0

∂ − ∂x

∂ ∂x

0



 ,

V=

{0}3×3

{0}3×3

ε−1 {R}3×3

{0}3×3

 (2b)

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(2c)

where ε,  are the dielectric constant and magnetic permeability of medium, respectively.

H.W. Yang, H. Song / Optik 124 (2013) 1716–1720

1717

2.1. Discretization in the time direction n+s/5

The time evolution of (2a) from t = 0 to t = t can be written as

  E

H

Hy

n+(s−1)/5

(i, j) = Hy

  E

(t) = exp(t(U + V ))

H

(0)

exp(t(U + V )) =

t 1 c √ + r l 0 ε0

(3)

where exp(t(U + V)) is the time evolution matrix of Maxwell’s equations. With symplectic propagator technique, there have m stage p-order approximation of exp(t(U + V)), that is: m

exp(dl tV ) exp(cl tU) + O(t

p+1

)

(i + 1/2, j)

n+s/5

Hx

n+(s−1)/5

(i, j) = Hx





∂Ez n+(s−1)/5 ∂x

(i + 1/2, j)

(i, j + 1/2)

t 1 c √ + r l 0 ε0







∂Ez n+(s−1)/5 ∂y

(8b)

 (i, j + 1/2)

(8c)

(4)

l=1

where cl and dl is the symplectic operator or the symplectic propagator coefficients. That can prove it. When Ul = Vl = 0 (l ≥ 2), using Taylor expansion technique (4) can be written as: exp(t(U + V )) =

m

(I + dl tV )(I + cl tU) + O(t p+1 )

2.2. Discretization in the space direction A function of space and time evaluated the true solution F at an approximate discrete point (ix, jy, kz) in the nth time step, can be notated as Fn (i, j, k) = F(ix, jy, kz ; nt). For the spatial direction, the high-order accuracy discrete scheme is:



q/2−1

=

Wr

r=0



(h)

F n+s/m (h − r − 1/2) − F n+s/m (h + r + 1/2) ı

+ O(ıq )

(6)

where ı = x,y, h = i,j and Wr are the coefficients of spatial differences. We note that at the m stage in time domain and the q-order in spatial domain of SFDTD method is SFDTD (m,q). At the same time, based on the two-dimensional TM case in ideal medium, Maxwell’s equations can be written as: ∂Hx 1 ∂Ez =−  ∂y ∂t

(7a)

∂Hy 1 ∂Ez =  ∂y ∂t

(7b)

1 ∂Ez = ε ∂t

∂Hy ∂Hx − ∂x ∂y

n+s/5

(9a)

ω(ω + ie )

1−



2 ωpm

(9b)

ω(ω + im )

where  e and  m are the loss coefficient. In order to explore the positive and negative of ε,  relationship with the medium of electromagnetic properties, it is discussed of ε and  role in the Maxwell equation. According to (1), introduction B(ω) = (ω)H(ω), D(ω) = ε(ω)E(ω), then use of the Z changes, it will transform the dielectric constant and magnetic permeability rate of a medium with frequency-dependent from frequency domain to the Z domain, the time domain convolution integral into a multiplicative of Z domain. Some modifications are made for formula (7a)–(7c). We select 2-order in time domain and 4-order in spatial domain of the 2 stage SFDTD (2,4) simulation, that can be further calculated the SFDTD differential equations of the Left-handed materials medium. The following only shows the electric field SFDTD difference equation [11–13], the magnetic field can be deduced according to this method. n+s/2

Dz

n+(s−1)/2

(i, j) = Dz +

1 t d √ εr l 0 ε0

−

Ez

(i, j)

∂Hx n+(s−1)/2 ∂y

n+(s−1)/2

(i, j) = Dz

·Ez



(i, j) −

∂Hy n+(s−1)/2 ∂x



 (i + 1/2, j)



(i, j + 1/2)

(10a)

2 · t ωpe

z −1 (ee t − 1)

e

(1 − z −1 )(1 − ee t z −1 )

(i, j)

(10b)

(7c)

(i, j)

1 t + dl √ εr 0 ε0 −

1−

n+(s−1)/2



n+(s−1)/5

(i, j) = Ez



2 ωpe

(ω) = 0

n+s/2

As follows, here are the (s − 1)th stage to sth stage discrete iteration formats of SFDTD method. Ez

ε(ω) = ε0

As for the coefficient of the specific values, it refers to the literature. The constant coefficient of the propagator obeys symmetry relations cl = cm−l+1 (1 ≤ l ≤ m), dl = dm−l (1 ≤ l ≤ m − 1) and dm = 0.

∂F n+s/m ∂ı

Drude model is chosen as a model of electromagnetic parameters of left-handed materials, relative permittivity and relative permeability model with Drude model [10]:

(5)

l=1

3. Left-handed material model



∂Hx n+(s−1)/5 ∂y

∂Hy n+(s−1)/5 ∂x



 (i + 1/2, j)



(i, j + 1/2)

(8a)

4. SFDTD advantage validation In vacuum, a case of one-dimensional Gaussian pulse Ei (t) = exp[− 0.5(t − t0 )2 / 2 ] is chosen to simulate the propagation process. The space step length is 1 cm, 10 grid is divided per wavelength, CFL = 0.5. Fig. 1(a) shows the analytical solution of the Ei (t) and FDTD method, SFDTD (2,4) method comparison chart when the time step in the 400th. Fig. 1(b) shows the analytical solution of the Ei (t) and FDTD method, SFDTD (2,4) method comparison chart when the time step in the 3000th.

1718

H.W. Yang, H. Song / Optik 124 (2013) 1716–1720

Fig. 2. Refraction models of different materials and spatial modeling. Fig. 1. FDTD method and SFDTD (2,4) method with the analytical solution of the electric field curve comparison.

It can be seen early in the calculation; different algorithms on the solution of pulse, waveforms are in good inosculation. But as time step increases, the waveform of FDTD method distorts beginning, and has the energy dissipation. The amplitude decreases. In the simulation of pulse propagation, whether pre-or post-calculation, the analytical solution of the pulse inosculate very well, and the amplitude of electric field is no obvious decay, the peak is always in one nearby. Although SFDTD (2,4) method has slightly energy loss, but the basic structure and the peak also remain at one near. It explains symplectic method maintaining good conservation within the Hamilton system. The energy conservation in wave propagation proved the stability and high accuracy of symplectic method. 5. The simulation results and analysis of left-handed materials electromagnetic characteristics Negative refraction analysis of the left-handed materials, and modeling in two-dimensional space is shown in Fig. 3. Drude model structure is selected by left-handed materials. SFDTD (2,4) simulates calculation region for 220 × 430. Left-handed materials parameters are = 0.01 m, f = 30 GHz, x = y = /100, √ ωp = 2 f. When refractive index is n = −1, that is ωpe = ωpm = 2ωp , loss coefficient  = 3.75 × 10−4 ωp ,  e =  m =  . Around the space, sets of 10 grids of perfectly matched layer (PML) as absorbing boundary. Left-handed material thickness is d2, added source distance from the front of the left-handed materials to d1 (Fig. 2). 5.1. Negative refraction of left-handed materials validation This section, discusses the negative refraction of Gaussian beam from free space oblique incidence to the left-handed material slab. When refractive index is n = −1, selects at the 100th mesh add the left-handed medium plate. The center of wave source from the left-hand media interface is the 80 mesh, that is d1 = 8 mm. Thickness left-handed media is 100 mesh, that is d2 = 10 mm. By the refractive index n = −1, the left-handed material layer and the free space l to match. For any angle of incidence cases, the reflection

and transmission coefficients of incident wave were R = 0, T = 1. According to Snell theorem, refraction angle can be drawn:  t = −  i . This section selects the incidence angle:  i = arctan(1/3). The results shown in Fig. 3. Fig. 3(a) for the incident wave electric field distribution of the SFDTD (2,4) method in vacuum. Fig. 3(b) for the refraction of electric field distribution of the SFDTD (2,4) method after the incident wave through left-handed materials. Fig. 3 compares 3050th time step, the Gaussian beam in free space and left-handed materials layer, the transmission electric field distribution. Refracted beam with the incident beam are located in the same side of normal in Fig. 3(b). They satisfy Snell’s law, that shows the Gaussian beam in left-handed materials and air occurs negative refraction at the interface. Among them, the continuous Gaussian beam oblique incidence to the left-handed material slab, the divergence of the beam is changed. 5.2. Left-handed material flat collection effect verification According to left-handed materials with negative refractive index properties, left-handed material slab has a collection effect. When a bunch of divergence wave from the ordinary medium plate by left-handed materials, this occurs the collection effect. And when the left-handed dielectric layer thick enough, it will form a convergence point within the left-handed dielectric layer. That is called confocal phenomenon. Suppose wave source normal incidence from the boundary of the total-scattering field. When the other initial conditions is as the same as previous section, two cases that the left-handed materials are 100 meshes and 180 meshes thickness are simulated at Fig. 4. Fig. 4(a) shows when t = 3050th step the wave transmission in free space. Comparison of Fig. 4(b) and (c), when the wave propagation in the left-handed dielectric layer, the wave reaches a steady state, the beam shape becomes the convergence shape from to diverging gradually. 5.3. Validation phase compensation effect of left-handed materials Fig. 5 shows phase compensation model for the left-handed materials. Phase compensation effect of the left-handed materials

H.W. Yang, H. Song / Optik 124 (2013) 1716–1720

1719

Fig. 3. Negative refraction of SFDTD (2,4) oblique incident Gaussian beam in lefthanded materials and free space.

is another form for its backward wave property. Electromagnetic wave propagation in left-handed materials caused by the phase advance can be used to make up the phase lag for dissemination in the general media generated. In Fig. 6, sets slab thickness of the left-handed materials and the right hand materials both account for 50 meshes size. Refractive index of left-handed materials is taken as −2 and −3, the refractive index of the righthanded material corresponding for the 2 and 3. Wave source is normal incidence from the boundary of the total-scattering field, the other initial conditions as the same as previous section. By the equation p = |n1 | k0 d1 − |n2 | k0 d2 , phase variation is zero. That is a phase compensation effect. When the time step for the 3050th, Fig. 6 shows the electric field intensity distribution with Gaussian beam through different refractive index of the left and right materials composed of phase compensation structure. Fig. 6 shows that after the Gaussian beam through the right-hand material flat, it will be broadening with the increase of propagation distance. However, after the beam through the left-handed material slab, it is re-assembled again. Finally, it spreads with the original trend in the free space. In addition, the left and the right materials with different refractive index, and between which wave propagation time is difference. Therefore,

Fig. 4. Negative refraction effect of SFDTD (2,4) method Gaussian beams at different thickness of the left-handed materials.

you can change the refractive index of the left-handed materials and the right-hand materials to control the propagation time of the Gaussian beam through the phase compensation structure slab.

1720

H.W. Yang, H. Song / Optik 124 (2013) 1716–1720

6. Conclusion

Fig. 5. Phase compensation structures of left-handed materials and right-handed materials.

This is the first use of symplectic finite-difference time-domain (SFDTD) (2,4) method of the Drude model to simulate the anomalous electromagnetic properties in left-handed materials. The method is based on dynamic system expressed in Hamilton, it can maintain the energy conservation (i.e. Hamilton volume conservation), and the corresponding mathematical formulation for the evolution of the system state is symplectic transformation. Thus in the calculation, the system state of the discretization equations should also be symplectic transformation. Because the system with symplectic algorithm can guarantee the time evolution is always a symplectic transformation, always keep the basic features of Hamilton system, which ensures that symmetry and conservation of the numerical method. The symplectic finite-difference time-domain (SFDTD) (2,4) method were compared with the general FDTD method, the results show that SFDTD (2,4) has the advantage of stability and high precision. This paper uses SFDTD (2,4) theory to analyze propagation characteristics of the electromagnetic wave through the left-handed materials, which proves the negative refraction effect, collection effect and phase compensation effect of the left-handed materials. References

Fig. 6. Phase compensation effects of the different refractive index plate with SFDTD (2,4) method.

[1] V.G. Veselago, The electrodynamics of substances with simultaneously negative values of ε and , Sov. Phys. Usp. 10 (January (4)) (1968) 509–514. [2] I.V. Lindell, S.A. Tretyakov, K.I. Nikoskinen, S. Ilvonen, BW media-media with negative parameters, capable of supporting backward waves, Microw. Opt. Technol. Lett. 31 (October (2)) (2001) 129–133. [3] M. Kusaf, A.Y. Oztoprak, D.S. Daoud, Optimized exponential operator coefficients for symplectic FDTD method, IEEE Microw. Wireless Comp. Lett. 15 (February (2)) (2005) 86–88. [4] R. Rieben, D. White, G. Rodrigue, High-order symplectic integration methods for finite element solutions to time dependent Maxwell equations, IEEE Trans. Antennas Propag. 52 (August (8)) (2004) 2190–2195. [5] I. Saitoh, Y. Suzuki, N. Takahashi, The symplectic finite difference time domain method, IEEE Trans. Magn. 37 (September (5)) (2001) 3251–3254. [6] D.B. Ge, Y.B. Yan, Finite-difference Time-domain Method for Electromagnetic Waves, 2nd edition, Xidian University Press, Xi’an, 2005. [7] D.R. Smith, N. Kroll, Negative refractive index in left-handed materials, Phys. Rev. Lett. 85 (October (14)) (2000) 2933–2936. [8] H. Chen, L. Ran, J. Huangfu, X. Zhang, K. Chen, Left-handed materials composed of only s-shaped resonators, Phys. Rev. E 70 (November (5)) (2004), pp. 0576051–057605-4. [9] D.R. Smith, W.J. Padilla, D.C. Vier, S.C. Nemat-Nasser, S. Schultz, Composite medium with simultaneously negative permeability and permittivity, Phys. Rev. Lett. 84 (May (18)) (2000) 4184–4187. [10] L.-l. Jiang, J.-f. Mao, X.l. Wu, Symplectic finite-difference time-domain method for Maxwell equations, IEEE Trans. Magn. 42 (August (8)) (2006) 1991–1995. [11] M. Kusaf, A.Y. Oztoprak, Higher stability limits for the symplectic FDTD method by making use of Chebyshev polynomials, IEEE Microw. Wireless Comp. Lett. 16 (November (11)) (2006) 579–581. [12] T. Hirono, W. Lui, S. Seki, Y. Yoshikuni, A three-dimensional fourthorder finite-difference time-domain scheme using a symplectic integrator propagator, IEEE Trans. Microw. Theory Tech. 49 (September (9)) (2001) 1640–1648. [13] W. Sha, Z.X. Huang, M.S. Chen, X.L. Wu, Survey on symplectic finite difference time-domain schemes for Maxwell’s equations, IEEE Trans. Antennas Propag. 56 (February (2)) (2008) 493–500.