High-order unified symplectic FDTD scheme for the metamaterials

Computer Physics Communications 183 (2012) 1192–1200

Contents lists available at SciVerse ScienceDirect

Computer Physics Communications www.elsevier.com/locate/cpc

High-order unified symplectic FDTD scheme for the metamaterials Xingang Ren a , Zhixiang Huang a,∗ , Xianliang Wu a,b , Silong Lu a , Hui Wang a , Lei Wu a , Shen Li a a b

Key Lab of Intelligent Computing and Signal Processing, Anhui University, Ministry of Education, Hefei, China School of Electronic and Information Engineering, Hefei Normal University, Hefei, China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 13 August 2011 Received in revised form 16 January 2012 Accepted 29 January 2012 Available online 1 February 2012 Keywords: Left-handed materials (LHMs) Unified symplectic finite-difference time-domain (US-FDTD) Split perfectly matched layers (SPML) Metamaterials (MTMs)

A high-order unified symplectic finite-difference time-domain (US-FDTD) method, which is energy conserved, for modeling the metamaterials is proposed. The lossless Drude dispersive model is taken into account in US-FDTD scheme, and the detailed formulations of the proposed US-FDTD method are also provided. The high-order split perfectly matched layers (SPML) are used as the absorbing boundary conditions (ABCs) to terminate the computational domain. The analysis of Courant stability and numerical dispersion demonstrate that US-FDTD scheme is more efficient than the traditional time domain numerical methods. Focusing and refocusing of the electromagnetic wave in target detection is validated using the normal incident Gaussian beam with a matched slab. Oblique incidence results associated with the inverse Snell effect and the phase compensation effect of the composite slab further demonstrated the efficiency of the method. Numerical results for a more realistic structure are also included. All the results agree well with the theoretical prediction. The method proposed here can be directly put into using as a time-domain full-wave simulation tool for applications in metamaterials. © 2012 Elsevier B.V. All rights reserved.

1. Introduction In 1968, Veselago predicated a kind of composited material named left-handed materials (LHMs), metamaterials (MTMs), or double-negative-index materials (DNGs), of which the permittivity ε and the permeability μ are both negative. It was characterized by some rather unusual properties, such as inverse Snell effect, an inverse Doppler shift, and backward-directed Cherenkov radiation [1]. In 2000, Pendry et al. [2] demonstrated that materials with an array of split ring resonators (SRRs) produce negative permeability over certain frequency bands. In 2000, Smith et al. [3] demonstrated the experimental existence of LHMs for the first time, combining an array of SRRs interspersed with a 2D array of wires. And also in 2000, Pendry [2] demonstrated that LHMs focus light “perfectly” even when they are in the form of a parallel-sided slab, which is unconventional to a classic lens. As it is known that metamaterials have the potential to control and manipulate the wave propagation in a manner that differs from the conventional materials because of their periodic nature. The impact of metamaterials could be enormous. One of the usages may enable us to create highly directional antennas, enhance the performance of small antennas and design highly integrated transceiver systems that can be packaged in a limited space. Although the practical applications of the metamaterials are currently limited by their loss characteristics and operational bandwidths, it is very mean-

*

Corresponding author. E-mail address: [email protected] (Z. Huang).

0010-4655/$ – see front matter doi:10.1016/j.cpc.2012.01.021

© 2012

Elsevier B.V. All rights reserved.

ingful to develop efficient modeling tools to quantify and represent the characteristics of the metamaterials and their behavior, especially when they form an integral part of complex electromagnetic system. Finite-difference time-domain (FDTD) is a timedomain technique, which provides transient electromagnetic field directly with a flexible choice of excitations and yields wideband result generated by a single simulation, and it is also a grid-based on numerical technique, which enables us to model objects with arbitrary material fillings, the codes can be easily parallelized. Ziolkowski [4] used the dispersive FDTD technique [5] as a tool to study the properties of metamaterials both analytically and numerically. Nowadays there are mainly three methods to obtain the numerical solutions of metamaterials. The auxiliary differential equation (ADE) [6] technique transforms the frequency-dependent constitutive relation in the time domain, by inverse Fourier transform, obtaining an ordinary differential equation. The Z-transform method [7] concludes in a similar differential equation, assuming the complex permittivity in the Z-domain to be a transfer function. The convolution integral of the recursive convolution (RC) [8] formulation corresponding to the time-domain constitutive relation is approximated by a discrete summation which is then properly calculated using a recursive procedure. The aforementioned dispersive FDTD methods approximate the time derivative with a low-order approximation (no more than second-order). And the conventional dispersive FDTD approach suffers from numerical inaccuracies when evanescent waves must be included in the modeling of the structure of interest. Furthermore, existing commercial software often does not offer the flexibility to modify in the algorithm used in their codes. So a high-order approximation in the

X. Ren et al. / Computer Physics Communications 183 (2012) 1192–1200

temporal domain is highly required. There were several methods proposed to overcome this shortcoming [9]. Fang [10] proposed the high-order FDTD(4, 4) method. Yet, the method is difficult to handle material interface for modeling the three-dimensional complex objects. The high-order Runge–Kutta (R–K) approach was introduced in [11,12] to improve the accuracy along the time direction. However the approach is dissipative and needs large amount memory. Another method is the alternating direction implicit FDTD (ADI-FDTD) algorithm [13–15]. Although it saves CPU time owing to unconditional stability, undesirable numerical precision and dispersion will happen once the high Courant–Friedrichs–Levy (CFL) number is adopted. Hence developing an energy-preserving highorder scheme both in temporal and spatial domain is of great necessity. The US-FDTD technique is introduced in this paper to model the metamaterials, based on the implementation of dispersive FDTD approach in conjunction with operator split technique. This operator split formulation represents the first of its kinds in the dispersive FDTD method. This technique can also be used as a tool to investigate the potentials application of general single and double dispersive materials, as well as efficiently explore the complicated and large scale electromagnetic systems. Especially, this technique can be reduced to the normal symplectic FDTD method that is flexible to deal with the composite materials. The basic idea of the US-FDTD is that Maxwell’s equations can be rewritten as an infinite-dimensional Hamiltonian system [16]. As is known, the Hamiltonian system preserves the energy and has an advantage of stability and low numerical dispersion in the numerical computation, especially under long term simulation. Since the proposed US-FDTD technique is a kind of the symplectic schemes, a stable and accurate solution can be obtained by using the symplectic schemes, which also preserves the energy of the electromagnetic system constant. The perfect matched layers (PML) [17] technique is very efficient for absorbing electromagnetic wave and solving the unbounded problems. The consistent high-order split perfect matched layers (SPML) [18] with a low reflectivity are used as the absorbing boundary conditions (ABCs) in the proposed US-FDTD scheme. The proposed US-FDTD can be used as a tool to gain a more efficient and accurate characterization of material properties. Here, in Section 2, the detail and theoretical formulae of the highorder unified symplectic finite-difference time-domain (US-FDTD) method are constructed, which have fourth-order accuracy both in temporal and spatial domain, to simulate the different problems in metamaterials. The lossless Drude dispersive model is introduced in US-FDTD method as an illustration. The detailed iterative formulations of the scheme are also presented. In Section 3, the stability and low numerical dispersion are compared with the other method firstly, then the numerical experiments of the perfect lens imaging system are also investigated, which shows that the USFDTD method results agree well with the theoretical analysis that Pendry predicted. Then the interesting effects of n = −1 for the normal and oblique incidence cases are also demonstrated. The phase compensation effect of the composite slab is also included in the last part. Although the numerical experiments are basic properties of the metamaterials which are idealized, the proposed US-FDTD can also be used to model the complicated structure and solve the realistic problem. 2. Formulations The lossless Drude model [4] is used to describe the metamaterials as:





ω2pe ω2   ω2pm μr (ω) = μr0 1 − 2 ω εr (ω) = εr0 1 −

(1a) (1b)

1193

where ω pe and ω pm are the plasma frequency of the material. Considering the constitutive relations then Maxwell’s equations in lossless and sourceless medium of the frequency domain can be written as:

˜ = −∇ × E˜ j ωμ0 μr (ω)H

(2a)

˜ j ωε0 εr (ω)E˜ = ∇ × H

(2b)

˜ H ˜ are the electric and magnetic field in the frequency where E, domain, ∇ is the curl operator, ε0 , μ0 are the vacuum permittivity ¯ 0 = μ0 μr0 for simplicity. and permeability, and let ε¯ 0 = ε0 εr0 , μ Substituting the Drude form of (1) into Eq. (2a) and (2b), we get

˜ =− j ωH j ωE˜ =

1

μ¯ 0

1

ε¯ 0

∇ × E˜ − K˜

(3a)

˜ − ˜J ∇ ×H

(3b)

˜ and ˜J are the auxiliary variables with Here K

ω2pm ˜ H ω ω2 ˜J = − j pe E˜ ω ˜ = −j K

(4a) (4b)

Applying the inverse Fourier transform, (3) and (4) can be rewritten in the time domain as

∂H 1 =− ∇ ×E−K ∂t μ¯ 0 ∂E 1 = ∇ ×H−J ∂t ε¯ 0 ∂K = ω2pm H ∂t ∂J = ω2pe E ∂t

(5a) (5b) (5c) (5d)

Eq. (5) can be expressed in a compact matrix form as

∂Ψ = AΨ ∂t

(6)

where Ψ = [H, E, K, J] T with the matrix

⎛ 0 − μ¯R0 ⎜ R 0 ε¯ 0 A=⎜ ⎝ ω2 I 0 pm ⎛

0

0 R = ⎝ ∂∂z

− ∂∂y

ω2pe I − ∂∂z 0 ∂ ∂x

−I

0 ⎞

0

−I ⎟ ⎟ 0 ⎠

0 0 ∂ ∂y − ∂∂x

(7)

0

⎞ ⎠

(8)

0

where I is the 3 × 3 identity matrix, 0 is the 3 × 3 null matrix and R = ∇× is the 3 × 3 matrix representing the curl operator. The time evolution of the electromagnetic field from t = 0 to t = t can be accurately obtained by using the exponential operator exp(At ) as

Ψ (t ) = exp(At )Ψ (0)

(9)

However due to the exponential operator exp(At ) cannot be evaluated at any t analytically. Fortunately, we can approximate exp(At ) to the pth order in the product form [19,20]:

1194

X. Ren et al. / Computer Physics Communications 183 (2012) 1192–1200

∂ f n+l/m

Table 1 Coefficients of the propagators. l

cl

dl

1 2 3 4 5

0.17399689146541 −0.1203850412143 0.89277629949778 −0.1203850412143 0.17399689146541

0.623379324513220 −0.12337932451322 −0.12337932451322 0.62337932451322





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

=

m

0



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



(10)

l =1

A=U+V ⎛0 − R μ¯ 0 ⎜0 0 U=⎝ 0 0 0 ω2pe I ⎛ 0 0 ⎜ ε¯R 0 0 V=⎜ ⎝ ω2 I 0 pm 0 0

(11)

−I 0 ⎞ 0 0⎟ ⎠

(12)

0 0 0 0 ⎞ 0 0 0 −I ⎟

0 0

⎟

(13)

exp(dl Ut) = I12×12 + dl Ut

(14a)

exp(cl Vtt ) = I12×12 + cl Vt

(14b)

Therefore, the symplectic scheme is applicable to Maxwell’s equations in temporal domain of the metamaterials. The detailed expressions of H and E at the l-stage are deduced as

Hn+l/m = Hn+(l−1)/m + cl t ×

−1  · En+(l−1)/m μ¯ 0

J

=J

, j, k   1 n+(l−1)/m i + , j, k = Ex 2     1 1 1 α y1 × H nz +l/m i + , j + , k + ε¯ 0 2 2   1 1 n+l/m i + , j − ,k − Hz 2 2    1 1 n+l/m i + , j, k + − αz1 × H y 2

(15a)

+ cl t ω

1

ε¯ 0

(15b)

 · Hn+l/m − Jn+l/m

(15c)

Kn+l/m = Kn+(l−1)/m + dl t ω2pm Hn+l/m

(15d)

where  is the discrete form of the curl operator R. Using the notation of Yee grid [20] for discretizing computational domain into grids and the space increments in the x-, y-, and z-directions are described as x,  y,  z. The function f at the l-stage of n time step can be expressed as





f (x, y , z; t ) = f i x, j  y , k z; (n + τl )t = f n+l/m i , j ,k

(16)

In this paper, the fourth-order finite difference, and the staggered Yee grid are used to approximate the first-order spatial derivatives. For example, ∂ f /∂ x can be approximated as

n+l/m ∂ f

∂x

2

 i+

1

1



2

2



n+l/m

2



1

3

n+l/m

27( f |

i + 12 , j ,k

n+l/m

− f|

i − 12 , j ,k

n+l/m

)−(f|

24x

i + 32 , j ,k

n+l/m

− f|

i − 32 , j ,k

) (17)



i + , j, k + − αz2 × H y 2 2   1 3 n+l/m i + , j, k − − Hy 2

n+l/m Jx

2





1

i + , j, k (18a) − dl t 2     1 1 n+l/m n+(l−1)/m i + , j, k = J x i + , j, k Jx 2 2   1 n+(l−1)/m i + , j, k (18b) + cl t ω2pe E x 2

t

n+l/m Hx

 i, j +

−1

t dl × 24 y −1 t αz2 = dl × 24 z

α y2 =

,

8 y t 9 αz1 = dl × , 8 z 1

1



,k + 2   1 1 n+(l−1)/m i, j + , k + = Hx 2

2

2

  1 n+(l−1)/m β y1 × E z i , j + 1, k + μ¯ 0 2   1 n+(l−1)/m i, j, k + − Ez 2    1 n+(l−1)/m i, j + , k + 1 − βz1 × E y

+

1





n+(l−1)/m

− Ey

 i, j +

1



2

,k 2   1 n+(l−1)/m i , j + 2, k + + β y2 × E z 2   1 n+(l−1)/m i , j − 1, k + − Ez 

i , j ,k

2

, j, k − 2    1 3 n+l/m i + , j + ,k + α y2 × H z 2 2   1 3 n+l/m i + , j − ,k − Hz

9

2 n+(l−1)/m pe E

En+l/m = En+(l−1)/m + dl t ×

≈

i+



1

α y1 = dl ×

− Kn+(l−1)/m n+(l−1)/m



n+l/m

where U and V are noncommuting operators, and cl , dl are the symplectic integrators which values can be found in Refs. [20,21], here we used the coefficient as shown in Table 1. Thus the 5-stage and fourth-order symplectic scheme of the lossless Drude model is constructed. In view of Uα = Vα = 0 (α  2), the exponential operators exp(cl Ut ) and exp(dl Vt ) can be computed analytically by Taylor series expansion as

n+l/m

n+l/m Ex

− Hy

0 ⎠ 0

∂ f n+l/m

Similarly, ∂ y |i , j ,k and ∂ z |i , j ,k can be constructed, respectively. The detailed expressions of the E x , H x , J x , and K x components at the l-stage in the US-FDTD scheme, without the absorbing boundary conditions, are as follows:

2

(19)

X. Ren et al. / Computer Physics Communications 183 (2012) 1192–1200

1195

Table 2 Stability number of different time-domain methods.

 − βz2 ×

Methods

FDTD(2, 2)

FDTD(2, 4)

US-FDTD(4, 2)

US-FDTD(4, 4)/US-FDTDC

2D 3D

0.707 0.577

0.606 0.495

1.062 0.867

0.910 0.743

n+(l−1)/m Ey

n+(l−1)/m

− Ey

 i, j +

 i, j + 1 2



1

,k + 2 

2

,k − 1

  1 1 n+l/m − cl t K x i, j + , k + 2 2   1 1 n+l/m i, j + , k + Kx 2 2   1 1 n+(l−1)/m i, j + , k + = Kx 2 2   1 1 n+l/m i, j + , k + + dl t ω2pm H x 2

9

t , 8 y t 9 βz1 = cl × , 8 z β y1 = cl ×

−1

2

t 24 y −1 t βz2 = cl × 24 z β y2 =

(20a)

(20b)

cl ×

(21)

3. Numerical validations

Fig. 1. (Color online.) Numerical dispersive curves for different time-domain methods.

For simplicity, only the two-dimensional (2D) TM polarization case is selected as the demonstration of the proposed US-FDTD scheme. However, our scheme can be easily extended to the general three-dimensional problems. 3.1. The advantage of the dispersive US-FDTD method The (D, E, B, H) and (E, J, H, M) schemes which are based on the ADE-FDTD method [22] both need to store the previous electric field and magnetic field value at every time step. The computational intensity will become especially larger when a large scaled electromagnetic problem is considered, while the dispersive USFDTD scheme may become more promising due to the low memory required. Using the same procedure as in Ref. [21], numerical stability of US-FDTD scheme can be obtained. The Courant stability number in the two-dimensional US-FDTD scheme is 0.910 compared with 0.707 for the traditional FDTD method. The number would be reduced to 0.606, if we apply fourth-order approximation to the spatial derivative with the traditional FDTD method, which is even lower than the US-FDTD method. Table 2 shows the Courant stability numbers of the common used FDTD and US-FDTD methods, and the first and second integer denote the accuracy in the time and space direction, respectively. Numerical dispersion of US-FDTD is compared with original FDTD, ADI-FDTD [23,24] and compact symplectic FDTD (CS-FDTD) [25] algorithms. First, the relative phase velocity error as a function of point-per-wavelength (PPW) for a plane wave traveling at θ = 60◦ and φ = 30◦ is shown in Fig. 1. Here the CFL number is set to be 0.4 which is at a low level. One can find out that CS-FDTD has the best numerical dispersion, and our scheme is superior to the original FDTD and ADI-FDTD approach. For a better knowledge of dispersion, the relative phase velocity error is firstly taken as a function of point-per-wavelength PPW and propagation angle φ , then as a function of point-per-wavelength PPW and CFL number in Fig. 2 and Fig. 3, respectively. At last, we show the relative phase velocity error as a function of the propagation angles θ and φ with

Fig. 2. (Color online.) Numerical dispersive curves for a plane wave traveling at θ = 60◦ , CFL = 0.5.

a fixed PPW and CFL number as shown in Fig. 4. From Figs. 1–4, the results indicate that the CS-FDTD has the best numerical dispersion, which may be due to the novel compact high-order spatial derivation. Our US-FDTD scheme is inferior to CS-FDTD just at a level of 5 dB but is superior to ADI-FDTD at a level of 20 dB. Considering the CS-FDTD is implicit scheme which costs huge computer resources due to matrix storage and computation in every time step while our scheme is explicit which saves the storage by substituting the previous electromagnetic field. So our scheme may

1196

X. Ren et al. / Computer Physics Communications 183 (2012) 1192–1200

in which T p is the period of one cycle, g on and g off are the window functions which are given by

g on (t ) = 10x3on − 15x4on + 6x5on g off (t ) = 1 −



10x3off

− 15x4off

+ 6x5off

(23a)



(23b)

where xon (t ) = 1 − (mT p − t )/mT p , xoff (t ) = [t − (m + n) T p ]/mT p . These smooth excitation functions generate minimal noise as the waves are introduced into the simulation region. The function g on (t ) goes smoothly from 0 to 1 in m-periods; the function g off (t ) goes smoothly from 1 to 0 in m-periods; so the signal maintains constant amplitude for its envelope over n-periods. For the cases considered below, we set the parameters m = 5 and n = 100. 3.3. Split perfect matched layer (SPML) The split perfect matched layer (SPML) is very efficient for absorbing electromagnetic waves and solving the unbound problems. The formulation of US-FDTD in lossy medium can also be derived in a similar manner. Fig. 3. (Color online.) Numerical dispersive curves for a plane wave traveling at θ = 60◦ , φ = 30◦ .

Hn+l/m = exp(−χ )Hn+(l−1)/m +

 ×

χ = cl ·

n+(l−1)/m

− · E

μ¯ 0

cl t (1 − exp(−χ ))

χ

− Kn+(l−1)/m

 (24a)

σ m t μ¯ 0

(24b)

Jn+l/m = Jn+(l−1)/m + cl t ω2pe En+(l−1)/m

(24c)

dl t (1 − exp(−ξ ))

En+l/m = exp(−ξ )En+(l−1)/m +

ξ 



 · Hn+(l−1)/m − Jn+(l−1)/m × ε¯ 0 σ t ξ = dl · ε¯ 0

(24d) (24e)

Kn+l/m = Kn+(l−1)/m + dl t ω2pm Hn+l/m

(24f)

where σ and σ m are electric and magnetic conductivity in the lossy region. So the formulae of the electric field component E xy and E xz in the SPML region [18,27] are described as following: Fig. 4. (Color online.) Numerical dispersive curves for a plane wave traveling at PPW = 10 and CFL = 0.4.

n+l/m

 i+

E xy

     1 −dl σz t n+(l−1)/m , j , k = exp E xy i + , j, k 2 ε¯ 0 2 1

be a more powerful and potential technique when manipulating the dispersive materials. In addition, we should mention that the proposed dispersive US-FDTD method is very flexible. We can simulate single dispersive medium by setting ω pe = 0 or ω pm = 0, and the method will become the normal symplectic FDTD method if ω pe = ω pm = 0 is adopted.

+

−dl σz t ε¯ 0 ) dl σz t ε¯ 0

1 − exp(

    1 1 n+l/m × α y1 × H z i + , j + ,k n+l/m

− Hz

i+

 + α y2 ×

3.2. Multiple cycle m–n–m pulses

n+l/m

The input time signals are multiple cycle m–n–m pulses [26] given by the expression:

⎧ g on (t ) sin(ω0 t ) 0  t < mT p ⎪ ⎪ ⎨ sin(ω0 t ) mT p  t < (m + n) T p (22) S (t ) = ⎪ g off (t ) sin(ω0 t ) (m + n) T p  t < (m + n + m) T p ⎪ ⎩ 0 (m + n + m) T p  t

− Hz n+l/m

E xz

 i+

1 2





, j , k = exp +

2



1 2

1

, j − ,k

n+l/m Hz

 i+

1 2

2



2



i+

1 2

3

3

, j + ,k 2



, j − ,k



2



(25)



1 −dl σz t n+(l−1)/m E xz i + , j, k ε¯ 0 2 −dl σz t ε¯ 0 ) dl σz t ε¯ 0

1 − exp(



X. Ren et al. / Computer Physics Communications 183 (2012) 1192–1200

1197

Fig. 5. (Color online.) Focusing action of matched LHMs slab.





× −αz1 × n+l/m

 i+

− Hy

 − αz2 × n+l/m

− Hy

n+l/m Hy

1 2

i+

1 2

i+

, j, k −

n+l/m Hy





 i+

1 2

, j, k −

1

1 2

, j, k +



2

, j, k + 3



2

3

1



2



2 (26)

where σz = σz (i + 1/2, j , k) is the local electric conductivity at point (i + 1/2, j , k) in the PML region, while the polynomial conductivities are employed varying from zero at vacuum layer interface to the maximum at the outside of the PML layer. And the expressions of other subcomponents can be constructed similarly. 3.4. Structure of the LHMs slab In this section, the center frequency was chosen to be f 0 = 30 GHz. For the cases of cylindrical wave and the oblique incidence beam, the simulation domain occupies 300 × 300, and the point-per-wavelength (PPW) was chosen to be PPW = 10, i.e.,  = λ0 /10. While for the normal incidence beam and phase compensation effect cases, the computational region was 360 × 360 with PPW = 30. Ten layers of high-order split perfectly matched layers (SPML) [18] as expressions (25) are employed to absorb the outgoing wave. The width and length of the slab are d1 and L, respectively. The source-to-slab distance is d0 = d1 /2, thus the conditions specified by the analytical solutions in Ref. [4] inside and outside of the LHMs slab structure are nicely met. When excited the multiple cycle m–n–m pulses located in front side of the LHMs slab, it can be focused at a focal point inside the LHMs lens and then refocused at another focal point outside the LHMs lens, as shown in Fig. 5. 3.5. Numerical results of the impedance matched LHMs slab 3.5.1. Normal incidence results √ In this case the parameters were ω p = 2ω0 , thus εr (ω0 ) = μr (ω0 ) = −1. Then the index of refraction is n = −1 [28–31]. Thus the LHMs slab is impedance matched with the free space (n = 1). First, the conventional FDTD and the proposed dispersive US-FDTD are used to simulate this matched LHMs slab with a cylindrical wave source. The LHMs slab is located in the middle of the domain, with I 0 = [120, 180] and J 0 = [20, 280] cells, corresponding to d1 = 60, L = 260. The line source is located at [ I s , J s ] = [90, 150], corresponding to d0 = 30. The stability numbers of the FDTD and US-FDTD are 0.5 and 0.8, respectively. Fig. 6 shows the distribution of the electric field at t = 75T p . It can be seen that the source point S can be focused at point F 1 inside the LHMs slab and refocused at point F 2 outside the LHMs slab, respectively. Especially, the US-FDTD with a larger CFL(= 0.8) can also get reasonable results, while traditional FDTD method will be divergent as expected if the CFL = 0.8 is selected.

Fig. 6. (Color online.) E z -field distribution at t = 75T p for the slab with n = −1, the space occupies 300 × 300 cells, PPW = 10, (a) is the ADE method with CFL = 0.5; (b) is the proposed US-FDTD method with CFL = 0.8. Table 3 Comparisons of different time-domain methods. Results

FDTD(2, 2)

FDTD(2, 4)

US-FDTD(4, 2)

US-FDTD

Physical time CFL PPW Spatial step No. of step Total run time Average time Memory

2500 ps 0.5 10 0.001 m 1500 35.3 s 0.0235 s 3.0 M (byte)

2500 ps 0.5 10 0.001 m 1500 41.0 s 0.0273 s 3.1 M (byte)

2500 ps 0.7 5 0.002 m 536 14.5 s 0.0271 s 0.9 M (byte)

2500 ps 0.7 5 0.002 m 536 18.7 s 0.0349 s 2.0 M (byte)

Additional numerical tests were done for further investigating the advantage of the US-FDTD method. The parameters are the same as those in Fig. 6, except the CFL number and the value of the PPW. Since all the numerical results and plots are the same as in Fig. 6, we will only focus on the results of the computational cost, which are listed in Table 3. It should be noted that the ADE-FDTD results will be divergent at CFL = 0.7 or become wrong with PPW = 5. As can be inferred from Table 3, the US-FDTD with coarse discretized grid and large CFL number can still obtain satisfactory results, which in turn proves to be a promising method, with advantages of high accuracy, low computational resources and facility of large domain and long-time simulation. Now the simulations with the index of refraction n = −2 and n = 0 [32] are also done for further checking the validity of the US-FDTD. The simulated time is 938 time steps, i.e., 75T p . Fig. 7a shows the distribution of the electric field for n = −2. In contrast to the case of n = −1, there is little focusing but a phase delay can be detected. Fig. 7b shows the distribution for the case of n = 0. The cylindrical wave can be converted to a plane wave when it passes through the LHMs slab, while the amplitude is relative small compared to the incident wave. These results would be anticipated configuration in metamaterials. For further investigating the physical meaning of n = −1, a Gaussian beam is also introduced [26]. The beam with a width of 30 grids propagates along the axis direction. The beam propagating in the free space is shown in Fig. 8a. The distribution of the electric field at t = 1510t clearly shows that the Gauss beam is

1198

X. Ren et al. / Computer Physics Communications 183 (2012) 1192–1200

Fig. 7. (Color online.) E z -field distribution after 75 periods, the parameters are the same as in Fig. 6 except for the index of the refraction (a) n = −2; (b) n = 0.

diverging as it propagates, and the electric field E z is decreasing as demonstrated in Fig. 8c. The interaction of the beam with an LHMs slab is also shown in Fig. 8b. The width and length of the slab are d1 = 180 and L = 320. The simulated time is 1510t again. It clearly reveals that the LHMs slab turns the diverging beam towards the beam axis, hence, acts as if the slab focuses the beam, corresponding to the amplitude of the electric field increasing as in Fig. 8d. Also, we found that the distribution of the electric field was symmetric at the front interface and the back interface of the slab, respectively. 3.5.2. Oblique incidence results The cases considered here deal with the issue of whether an LHMs slab can provide a negative refraction or not with oblique incidence [26,33]. The beam considered here is dealing with the index of refraction n = −1, −2 and 2. In all cases the angle of incidence of the beam is π /4 and the simulation time is 75T p . As shown in Fig. 9, the negative angle of refraction is clearly seen at time step t = 75T p . Though the case is an oblique incidence, the beam did focus when through the LHMs slab, and the reflected beam is much smaller in contrast to the transmitted beam, because the index of refraction of LHMs slab is matched with the free space. The Snell effect can be detected in Fig. 10a when the beam crosses the slab with n = 2. On the other hand, as indicated in Fig. 10b, the inverse Snell effect is clearly seen. 3.5.3. Phase compensation effect of the composite slab The composite structure is composed of an RHMs and an LHMs slab, and the wave vector k and the Poynting vector S in the free space and composite slabs are also depicted as in Fig. 11 [34,35]. The index of the refraction is n R = 2 and n L = −2 of the RHMs and LHMs slab, respectively. It is well known that a beam will expand in the RHMs slab and will be refocused in the LHMs slab. Though there should be a loss in amplitude in this process, the phase of the beam at the back face of the composite slab should be the same as the value at the front face. The distribution of the electric field at time steps 1510t is shown in Fig. 12a. The expansion of the beam in the RHMs slab and the refocusing in the LHMs slab are clearly demonstrated. In Fig. 12b, the E z -field at the front face and the back face of the

Fig. 8. (Color online.) E z -field distribution for a Gauss beam at t = 1510t with the US-SFDTD. The space occupies 660 × 360 grids with PPW = 30 and CFL = 0.8. (a) Gauss beam propagating in the free space and (b) there is a slab with a width d1 = 180 and length L = 320 in the space. (c, d) E z -field distribution along the beam axis. It’s noted that the area among the two red lines was the location of the LHMs slab in (c) and (d).

composite slab are also presented. As it can be seen from the figure, the results agree with each other qualitatively. A large error is also introduced near the interface due to the discontinuity of the fields. So a more effective technique should be adopted here to improve the efficiency of the US-FDTD scheme, which is the topic in our future research. 3.5.4. A more realistic plasmonic metamaterial structure In order to systematically investigate US-FDTD method, we apply our method to do the simulations for a more realistic plasmonic metamaterial structure [36] as shown in Fig. 13. The metamaterials studied here consist of periodic arrays of asymmetric split-ring slits, which have been successfully applied to switching,

X. Ren et al. / Computer Physics Communications 183 (2012) 1192–1200

1199

Fig. 9. (Color online.) The inverse Snell effect of the index of refraction n = −1 at t = 75T p corresponding to 938t. Fig. 12. (Color online.) Phase compensation effect of the composite slab, the LHMs slab changed to the same thickness composite slab (a) E z -field distribution at t = 1510t. (b) The E z -field distribution of the front and back face of the composite slab.

Fig. 13. (Color online.) Schematic of the plasmonic metamaterial. The yellow and the violet are the gold (Drude model) and substrate (n = 1.48), respectively. The corresponding parameters are h = 230 nm, h s = 180 nm, w = 65 nm, ax = a y = 470 nm, t = g = 170 nm and all the other parameters used are the same as in Ref. [36].

Fig. 10. (Color online.) The Snell effect of the index of refraction n = 2 (a) and n =

−2 (b) at t = 75T p corresponding to 938t.

nonlinear, and sensor applications. The incident wave propagates along the z-direction and has the electric field polarization along x- or y-direction. In x- and y-direction, a square unit is used with the unit size of 540 nm. The SPML absorbing boundary conditions are employed in z-direction. Along the unit cell boundaries in x- and y-directions, periodic boundary conditions are enforced to simulate the infinite periodic structure. To simulate this structure, a sufficiently large distance between the structure and the SPML layers should be enforced to prevent these instabilities from the interfering with the evanescent wave. Calculated spectra of transmission T , reflection R, and absorption A for the metamaterial are shown in Fig. 14. The results agree quantitatively well with the measured data as shown in Ref. [36] (Fig. 1), which indicates that our method can be applied directly in the simulation of the metamaterials. 4. Conclusions

Fig. 11. (Color online.) The structure for the phase compensation effect. The direction of the wave vector k and Poynting vector S in the free space and the composite slab.

In this paper, a high-order US-FDTD method for simulating the metamaterials was proposed. The theory and corresponding

1200

X. Ren et al. / Computer Physics Communications 183 (2012) 1192–1200

[12]

[13]

[14]

Fig. 14. (Color online.) Calculated spectra of transmission T , reflection R, and absorption A for (a) x polarization and (b) y polarization.

formulas are also presented. The numerical results of the proposed approach reveal a high stability and low numerical dispersion. Then both normal and oblique incidence cases are considered with different refractive index. Both the cylindrical wave and the Gauss beam can be refocused if the index of the refraction of the LHMs slab is n = −1 which is impedance matched with the free space. Also the index of refraction of n = 0, −2 are considered. The new and interesting effects are all demonstrated using our proposed US-FDTD scheme. Although the method presented here is mainly dealing with the two-dimensional problems, it can be easily extended to the general three-dimensional problems. Thus, our method can be directly put into using as a time-domain fullwave simulation tool for metamaterials applications. Acknowledgements The authors are greatly appreciating the valuable comments of the anonymous reviewers. The authors gratefully acknowledge the support of the National Natural Science Foundation of China (Nos. 60931002, 61101064), Distinguished Natural Science Foundation (No. 1108085J01), Universities Natural Science Foundation of Anhui Province (Nos. KJ2011A002, KJ2011A242) and the 211 Project of Anhui University. References [1] V.G. Veselago, The electrodynamics of substances with simultaneously negative values of permittivity and permeability, Soviet Physics 10 (4) (1968) 509–514. [2] J.B. Pendry, Negative refraction makes a perfect lens, Physical Review Letters 85 (18) (2000) 3966–3969. [3] D.R. Smith, W.J. Padilla, D.C. Vier, et al., Composite medium with simultaneously negative permeability and permittivity, Physical Review Letters 84 (18) (2000) 4184. [4] R.W. Ziolkowski, E. Heyman, Wave propagation in media having negative permittivity and permeability, Physical Review E 64 (5) (2001) 056625. [5] C.L. Holloway, P.M. McKenna, R.A. Dalke, et al., Time-domain modeling, characterization, and measurements of anechoic and semi-anechoic electromagnetic test chambers, IEEE Transactions on Electromagnetic Compatibility 44 (1) (Feb. 2002) 102–118. [6] R.M. Joseph, S.C. Hagness, A. Taflove, Direct time integration of Maxwell’s equations in linear dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses, Optics Letters 16 (18) (1991) 1412– 1414. [7] D.M. Sullivan, Frequency-dependent FDTD methods using Z transforms, IEEE Transactions on Antennas and Propagation 40 (10) (1992) 1223–1230. [8] R. Luebbers, F.P. Hunsberger, K.S. Kunz, et al., A frequency-dependent finitedifference time-domain formulation for dispersive materials, IEEE Transactions on Electromagnetic Compatibility 32 (3) (1990) 222–227. [9] N.V. Kantartzis, T.D. Tsiboukis, Higher-Order FDTD Schemes for Waveguide and Antenna Structures, Morgan & Claypool Publishers, San Rafael, CA, USA, 2006. [10] J. Fang, Time domain finite difference computation for Maxwell’s equations, Ph.D. dissertation, Univ. California, Berkeley, CA, 1989. [11] J.L. Young, D. Gaitonde, J.J.S. Shang, Toward the construction of a fourth-order difference scheme for transient EM wave simulation: staggered grid approach,

[15] [16]

[17] [18]

[19] [20]

[21]

[22]

[23]

[24]

[25]

[26] [27]

[28]

[29]

[30]

[31]

[32]

[33] [34]

[35]

[36]

IEEE Transactions on Antennas and Propagation 45 (11) (Nov. 1997) 1573– 1580. J.S. Shang, High-order compact-difference schemes for time-dependent Maxwell equations, Journal of Computational Physics 153 (2) (Aug. 10, 1999) 312–333. T. Namiki, A new FDTD algorithm based on alternating-direction implicit method, IEEE Transactions on Microwave Theory and Techniques 47 (10) (1999) 2003–2007. Z. Fenghua, C. Zhizhang, Z. Jiazong, Toward the development of a threedimensional unconditionally stable finite-difference time-domain method, IEEE Transactions on Microwave Theory and Techniques 48 (9) (2000) 1550– 1558. E.L. Tan, D.Y. Heh, ADI-FDTD method with fourth order accuracy in time, IEEE Microwave and Wireless Components Letters 18 (5) (May 2008) 296–298. W.E.I. Sha, Wu Xianliang, Huang Zhixiang, Chen Mingsheng, in: Vitaliy Zhurbenko (Ed.), The High-Order Symplectic Finite-Difference Time-Domain Scheme: Passive Microwave Components and Antennas, Intech, 2010. J.P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics 114 (2) (Oct. 1994) 185–200. T. Hirono, W.W. Lui, S. Seki, Successful applications of PML-ABC to the symplectic FDTD scheme with 4th-order accuracy in time and space, Microwave Symposium Digest, IEEE MTT-S International 3 (1999) 1293–1296. H. Yoshida, Construction of higher order symplectic integrators, Physics Letters A 150 (5–7) (1990) 262–268. W.E.I. Sha, Z.X. Huang, M.S. Chen, et al., Survey on symplectic finite-difference time-domain schemes for Maxwell’s equations, IEEE Transactions on Antennas and Propagation 56 (2) (2008) 493–500. T. Hirono, L. Wayne, S. Seki, et al., A three-dimensional fourth-order finitedifference time-domain scheme using a symplectic integrator propagator, IEEE Transactions on Microwave Theory and Techniques 49 (9) (2001) 1640–1648. Y. Zhao, P. Belov, Y. Hao, Accurate modelling of left-handed metamaterials using a finite-difference time-domain method with spatial averaging at the boundaries, Journal of Optics. A, Pure and Applied Optics 9 (9) (Sep. 2007) S468–S475. M. Kusaf, A.Y. Oeztoprak, Dispersion analysis of the ADI-FDTD and SFDTD methods, Turkish Journal of Electrical Engineering and Computer Sciences 16 (3) (2008) 201–206. N.V. Kantartzis, D.L. Sounas, C.S. Antonopoulos, et al., A wideband ADI-FDTD algorithm for the design of double negative metamaterial-based waveguides and antenna substrates, IEEE Transactions on Magnetics 43 (4) (Apr. 2007) 1329– 1332. J.Y. Wang, P. Liu, Y.L. Long, A compact symplectic high-order scheme for timedomain Maxwell’s equations, IEEE Antennas and Wireless Propagation Letters 9 (2010) 371–374. R. Ziolkowski, Pulsed and CW Gaussian beam interactions with double negative metamaterial slabs, Optics Express 11 (7) (2003) 662–681. W.E.I. Sha, Z.X. Huang, X.H. Wu, et al., Application of the symplectic finitedifference time-domain scheme to electromagnetic simulation, Journal of Computational Physics 225 (1) (2007) 33–50. X.D. Chen, T.M. Grzegorczyk, B.I. Wu, et al., Robust method to retrieve the constitutive effective parameters of metamaterials, Physical Review E 70 (1) (Jul. 2004). L.X. Ran, J. Huangfu, H. Chen, X.M. Zhang, K.-S. Cheng, T.M. Grzegorczyk, J.A. Kong, Experimental study on several left-handed metamaterials, Progress in Electromagnetics Research 51 (2005) 249–279. F. Mesa, M.J. Freire, R. Marques, et al., Three-dimensional superresolution in metamaterial slab lenses: Experiment and theory, Physical Review B 72 (23) (Dec. 2005). P.A. Belov, Y. Hao, Subwavelength imaging at optical frequencies using a transmission device formed by a periodic layered metal-dielectric structure operating in the canalization regime, Physical Review B 73 (11) (Mar. 2006). Y.G.G. Wang, H. Wang, On the size of left-handed material lens for near-field target detection by focus scanning, Progress in Electromagnetics Research 87 (2008) 345–361. K. Aydin, I. Bulu, E. Ozbay, Focusing of electromagnetic waves by a left-handed metamaterial flat lens, Optics Express 13 (22) (Oct. 31, 2005) 8753–8759. H. Luo, W. Hu, Z. Ren, et al., Focusing and phase compensation of paraxial beams by a left-handed material slab, Optics Communications 266 (1) (2006) 327–331. C. Caloz, A. Sanada, T. Itoh, A novel composite right-/left-handed coupled-line directional coupler with arbitrary coupling level and broad bandwidth, IEEE Transactions on Microwave Theory and Techniques 52 (3) (Mar. 2004) 980– 992. K. Tanaka, E. Plum, J.Y. Ou, T. Uchino, N.I. Zheludev, Multi-fold enhancement of quantum dot luminescence in a plasmonic metamaterial, Physical Review Letters 105 (2010) 2274031-4.