# FDTD Modeling of a Cloak with a Nondiagonal Permittivity Tensor

Authors: Naoki Okada; James B. Cole
Journal: ISRN Optics
 (2012)

---

## Abstract

We demonstrate a finite-difference time-domain (FDTD) modeling of a cloak with a nondiagonal permittivity tensor. Numerical instability due to material anisotropies is avoided by mapping the eigenvalues of the material parameters to a dispersion model. Our approach is implemented for an
elliptic-cylindrical cloak in two dimensions. Numerical simulations demonstrated the stable calculation and cloaking performance of the elliptic-cylindrical cloak.


---

## Body

## 1. Introduction

An optical cloak enables objects to be concealed from electromagnetic detection. Pendry et al. developed a method to design cloaks via coordinate transformations [1]. The coordinate transformation is such that light is guided around the cloak region. Material parameters (permittivity and permeability) can be obtained in the transformed coordinate system and put into Maxwell’s equations. This approach enables one to design not only cloaks but also other metamaterials that can manipulate light flow. For example, concentrators [2], rotation coatings [3], polarization controllers [4–6], waveguides [7–11], wave shape conversion [12], object illusions [13–15], and optical black holes [16, 17] have been designed. However, not many metamaterials have been realized in the optical region [18–25], because material parameters given by coordinate transformations have complicated anisotropies.Numerical simulations are useful to analyze complicated metamaterial structures. In this paper, we present a finite-difference time-domain (FDTD) analysis of a cloak. The FDTD method has gained popularity for several reasons: it is easy to implement, it works in the time domain, and its arbitrary shapes can be calculated [26–29]. FDTD modelings of cloaks with a diagonal (uniaxial) permittivity tensor have been demonstrated [30–38], but a cloak with a nondiagonal permittivity tensor has never been calculated by the FDTD method. The diagonal case can be stably calculated by mapping material parameters having values less than one to a dispersion model [31]. However, we found that mapping the nondiagonal elements to dispersion models causes the computation to diverge.In this paper, we analyze the numerical stability for a cloak with a nondiagonal permittivity tensor and derive the FDTD formulation. We apply our method to simulate light propagation in the vicinity of an elliptic-cylindrical cloak. To the best of authors’ knowledge, this is the first time that a cloak with a nondiagonal anisotropy has been calculated using the FDTD method.


## 2. Numerical Stability for Nondiagonal Permittivity Tensor

In the stability analysis, we confirm that the FDTD method for a cloak with a diagonal permittivity tensor cannot directly be extended to the nondiagonal case. Under a coordinate transformation for a cloak [39], material parameters can be expressed as
(1)ɛij=μij=±ggij,
where ɛij is the relative permittivity, μij is the relative permeability, gij is the metric tensor, and g=detgij. Because ɛij,μij are constructed from the symmetric metric tensor gij, they are symmetric. Consequently, ɛij,μij have real eigenvalues with orthogonal eigenvectors and are thus diagonalizable. The eigenvalues, λ, of ɛij,μij for an eigenvector V are defined by(2)ɛijV=μijV=λV.
The phase velocity of light in a material is given by c=c0/λ (c0= vacuum light speed), and the Courant-Friedrichs-Lewy (CFL) stability limit becomes(3)Δt≤λhc0d,
where Δt is the time step, h is the grid spacing, and d= 1, 2, and 3 dimensions. Since the FDTD stability depends on the eigenvalues of ɛij and μij, to analyze nondiagonal cases, we must first find the eigenvalues and diagonalize ɛij and μij. After the diagonalization, the FDTD method for diagonal cases [31–38] can be applied. For diagonal ɛij and μij, elements having values less than one are replaced by dispersive quantities to avoid violating the causality and numerical stability [40–44].In summary, the FDTD modeling for nondiagonalɛij and μij requires three steps:(1)
find the eigenvalues and eigenvectors and diagonalize the material parameters,(2)
map the eigenvalues having values less than one to a dispersion model,(3)
solve Maxwell’s equations using the dispersive FDTD method.


## 3. FDTD Formulation of the Elliptic-Cylindrical Cloak

Two designs of elliptic-cylindrical cloaks have been proposed. One has diagonalɛij and μij in orthonormal elliptic-cylindrical coordinates [45, 46], and in the other ɛij and μij are nondiagonal in Cartesian coordinates [47, 48]. We derive a FDTD formulation for the latter in the transverse magnetic (TM) polarization.


### 3.1. Diagonalization

Figure1 shows an elliptic-cylindrical cloak in Figure 1(a) Cartesian coordinates and Figure 1(b) transformed coordinates. The inner axis a, the outer axis b, and the perpendicular axes ka and kb are depicted. The elliptic-cylindrical cloak is horizontal when k>1, and vertical when k<1. In the cloak region, ka≤x2+k2y2≤kb, the material parameters are expressed by(4)ɛij=μij=[ɛxxɛxy0ɛxyɛyy000ɛzz],
where
(5)ɛxx=rr-ka+k2a2R2-2kar3(r-ka)r5x2,ɛxy=k2a2R2-ka(1+k2)r3(r-ka)r5xy,ɛyy=rr-ka+k2a2R2-2k3ar3(r-ka)r5y2,(6)ɛzz=(bb-a)2r-kar,
where r=x2+k2y2 and R=x2+k4y2. From (4) to (6) we can obtain three eigenvalues
(7)λ1=α-1α+1,λ2=1λ1,λ3=ɛzz,
where
(8)α=1+4r5(r-ka)k2a2R2(x2+y2).
Since ɛij is symmetric, it is diagonalized by the eigenvalue matrix Λ and its orthogonal matrix P as follows:
(9)ɛij=PΛPT,
where
(10)Λ=[λ1000λ2000λ3],(11)P=[ɛxyβ(λ2-ɛyy)β0(ɛyy-λ2)βɛxyβ0001],
where β=ɛxy2+(λ2-ɛyy)2.Elliptic-cylindrical cloak. (a) Cartesian coordinates:a, b are inner and outer axes; ka, kb are the perpendicular axes. (b) Transformed coordinates.
(a)(b)


### 3.2. Mapping Eigenvalues to a Dispersion Model

From (7) and (8), λ1 and λ3 have values less than one in the cloak region (ka≤r≤kb). Thus, λ1,λ3 must be replaced by dispersive quantities by using (for example) the Drude model
(12)λi=ɛ∞i-ωpi2ω2-jωγi,(i=1,3),
where ω is the angular frequency, ɛ∞i is the infinite-frequency permittivity, ωpi is the plasma frequency, and γi is the collision frequency. For simplicity, we consider the lossless case, γi=0. Then the plasma frequencies are given by ωpi=ωɛ∞i-λi, where ɛ∞i=max(1,λi).


### 3.3. FDTD Discretization

Using the diagonalized material parameters and eigenvalues mapped to the Drude model, we derive an FDTD formulation to solve Maxwell’s equations,(13)∂D∂t=∇×H,-∂B∂t=∇×E,
where D is the electric flux density, H is the magnetic field, B is the magnetic flux density, and E is the electric field. In the TM polarization, electromagnetic fields reduce to three nonzero components Ex, Ey, and Hz (Dx, Dy, and Bz). The D- and B-update equations are obtained using Yee algorithm [26–29] as follows:
(14)Dxn+1=Dxn+ΔthdyHzn+1/2,Dyn+1=Dyn-ΔthdxHzn+1/2,(15)Bzn+3/2=Bzn+1/2-Δth(dxEyn+1-dyExn+1),
where we simply write Dx(t=nΔt)→Dxn (n= integer) and dx, dy are the spatial difference operators defined by
(16)dxf(x,y)=f(x+h2,y)-f(x-h2,y),dyf(x,y)=f(x,y+h2)-f(x,y-h2).To find theE-update equations, we consider the relation
(17)D=ɛ0ɛijE,
where ɛ0 is the vacuum permittivity. From (9), we obtain
(18)ɛ0E=(ɛij)-1D=PΛ-1PTD.
Substituting (10) in (18) and multiplying λ1λ2 by both sides, we obtain
(19)ɛ0λ1λ2Ex=(λ1t22+λ2t12)Dx+t1t2(λ1-λ2)Dy,(20)ɛ0λ1λ2Ey=(λ1t12+λ2t22)Dy+t1t2(λ1-λ2)Dx,
where t1=ɛxy/β and t2=(λ2-ɛyy)/β. Substituting the Drude model for λ1 as shown in (12) and using the inverse Fourier transformation rule, -ω2→∂2/∂t2, (19) becomes
(21)ɛ0λ2(ɛ∞1∂2∂t2+ωp12)Ex=[(ɛ∞1t22+λ2t12)∂2∂t2+ωp12t22]Dx+t1t2[(ɛ∞1-λ2)∂2∂t2+ωp12]Dy.
For the discretization, we use the central difference approximation and the central average operator,
(22)∂2∂t2Exn=Exn+1-2Exn+Exn-1Δt2,Exn=Exn+1+2Exn+Exn-14.
The central average operator improves the stability and accuracy [40, 49, 50]. Similarly, Dx and Dy are discretized, and we obtain the Ex-update equation
(23)Exn+1=-Exn-1+2a1-a1+Exn+1ɛ0λ2a1+×[Dynb1+(Dxn+1+Dxn-1)-2b1-Dxn+c1+(Dyn+1+Dyn-1)-2c1-Dyn],
where Dyn+1,Dyn,Dyn-1 must be spatially interpolated due to the staggered Yee cell [31], and
(24)ai±=ɛ∞iΔt2±ωpi24,bi±=ɛ∞it22+λ2t12Δt2±ωpi2t224,ci±=t1t2ɛ∞i-λ2Δt2±ωpi2t224.
Similarly, the Ey-update equation is obtained by exchanging t1↔t2, Ex↔Ey, and Dx↔Dy.To find theHz-update equation, we consider the relation
(25)Bz=μ0ɛzzHz=μ0(ɛ∞3-ωp32ω2)Hz.
Analogously to the Ex field, the Hz-update equation can be obtained in the form
(26)Hzn+3/2=-Hzn-1/2+2a3-a3+Hzn+1/2+Bzn+3/2-2Bzn+1/2+Bzn-1/2μ0Δt2a3+,
where μ0 is the vacuum permeability.In summary, the electromagnetic fields are iteratively updated in the following sequence:(1)
update the components ofDn+1 according to (14),(2)
update the components ofEn+1 according to the sample given in (23),(3)
update the components ofBn+3/2 according to (15),(4)
update the components ofHn+3/2 according to (26).


## 3.1. Diagonalization

Figure1 shows an elliptic-cylindrical cloak in Figure 1(a) Cartesian coordinates and Figure 1(b) transformed coordinates. The inner axis a, the outer axis b, and the perpendicular axes ka and kb are depicted. The elliptic-cylindrical cloak is horizontal when k>1, and vertical when k<1. In the cloak region, ka≤x2+k2y2≤kb, the material parameters are expressed by(4)ɛij=μij=[ɛxxɛxy0ɛxyɛyy000ɛzz],
where
(5)ɛxx=rr-ka+k2a2R2-2kar3(r-ka)r5x2,ɛxy=k2a2R2-ka(1+k2)r3(r-ka)r5xy,ɛyy=rr-ka+k2a2R2-2k3ar3(r-ka)r5y2,(6)ɛzz=(bb-a)2r-kar,
where r=x2+k2y2 and R=x2+k4y2. From (4) to (6) we can obtain three eigenvalues
(7)λ1=α-1α+1,λ2=1λ1,λ3=ɛzz,
where
(8)α=1+4r5(r-ka)k2a2R2(x2+y2).
Since ɛij is symmetric, it is diagonalized by the eigenvalue matrix Λ and its orthogonal matrix P as follows:
(9)ɛij=PΛPT,
where
(10)Λ=[λ1000λ2000λ3],(11)P=[ɛxyβ(λ2-ɛyy)β0(ɛyy-λ2)βɛxyβ0001],
where β=ɛxy2+(λ2-ɛyy)2.Elliptic-cylindrical cloak. (a) Cartesian coordinates:a, b are inner and outer axes; ka, kb are the perpendicular axes. (b) Transformed coordinates.
(a)(b)


## 3.2. Mapping Eigenvalues to a Dispersion Model

From (7) and (8), λ1 and λ3 have values less than one in the cloak region (ka≤r≤kb). Thus, λ1,λ3 must be replaced by dispersive quantities by using (for example) the Drude model
(12)λi=ɛ∞i-ωpi2ω2-jωγi,(i=1,3),
where ω is the angular frequency, ɛ∞i is the infinite-frequency permittivity, ωpi is the plasma frequency, and γi is the collision frequency. For simplicity, we consider the lossless case, γi=0. Then the plasma frequencies are given by ωpi=ωɛ∞i-λi, where ɛ∞i=max(1,λi).


## 3.3. FDTD Discretization

Using the diagonalized material parameters and eigenvalues mapped to the Drude model, we derive an FDTD formulation to solve Maxwell’s equations,(13)∂D∂t=∇×H,-∂B∂t=∇×E,
where D is the electric flux density, H is the magnetic field, B is the magnetic flux density, and E is the electric field. In the TM polarization, electromagnetic fields reduce to three nonzero components Ex, Ey, and Hz (Dx, Dy, and Bz). The D- and B-update equations are obtained using Yee algorithm [26–29] as follows:
(14)Dxn+1=Dxn+ΔthdyHzn+1/2,Dyn+1=Dyn-ΔthdxHzn+1/2,(15)Bzn+3/2=Bzn+1/2-Δth(dxEyn+1-dyExn+1),
where we simply write Dx(t=nΔt)→Dxn (n= integer) and dx, dy are the spatial difference operators defined by
(16)dxf(x,y)=f(x+h2,y)-f(x-h2,y),dyf(x,y)=f(x,y+h2)-f(x,y-h2).To find theE-update equations, we consider the relation
(17)D=ɛ0ɛijE,
where ɛ0 is the vacuum permittivity. From (9), we obtain
(18)ɛ0E=(ɛij)-1D=PΛ-1PTD.
Substituting (10) in (18) and multiplying λ1λ2 by both sides, we obtain
(19)ɛ0λ1λ2Ex=(λ1t22+λ2t12)Dx+t1t2(λ1-λ2)Dy,(20)ɛ0λ1λ2Ey=(λ1t12+λ2t22)Dy+t1t2(λ1-λ2)Dx,
where t1=ɛxy/β and t2=(λ2-ɛyy)/β. Substituting the Drude model for λ1 as shown in (12) and using the inverse Fourier transformation rule, -ω2→∂2/∂t2, (19) becomes
(21)ɛ0λ2(ɛ∞1∂2∂t2+ωp12)Ex=[(ɛ∞1t22+λ2t12)∂2∂t2+ωp12t22]Dx+t1t2[(ɛ∞1-λ2)∂2∂t2+ωp12]Dy.
For the discretization, we use the central difference approximation and the central average operator,
(22)∂2∂t2Exn=Exn+1-2Exn+Exn-1Δt2,Exn=Exn+1+2Exn+Exn-14.
The central average operator improves the stability and accuracy [40, 49, 50]. Similarly, Dx and Dy are discretized, and we obtain the Ex-update equation
(23)Exn+1=-Exn-1+2a1-a1+Exn+1ɛ0λ2a1+×[Dynb1+(Dxn+1+Dxn-1)-2b1-Dxn+c1+(Dyn+1+Dyn-1)-2c1-Dyn],
where Dyn+1,Dyn,Dyn-1 must be spatially interpolated due to the staggered Yee cell [31], and
(24)ai±=ɛ∞iΔt2±ωpi24,bi±=ɛ∞it22+λ2t12Δt2±ωpi2t224,ci±=t1t2ɛ∞i-λ2Δt2±ωpi2t224.
Similarly, the Ey-update equation is obtained by exchanging t1↔t2, Ex↔Ey, and Dx↔Dy.To find theHz-update equation, we consider the relation
(25)Bz=μ0ɛzzHz=μ0(ɛ∞3-ωp32ω2)Hz.
Analogously to the Ex field, the Hz-update equation can be obtained in the form
(26)Hzn+3/2=-Hzn-1/2+2a3-a3+Hzn+1/2+Bzn+3/2-2Bzn+1/2+Bzn-1/2μ0Δt2a3+,
where μ0 is the vacuum permeability.In summary, the electromagnetic fields are iteratively updated in the following sequence:(1)
update the components ofDn+1 according to (14),(2)
update the components ofEn+1 according to the sample given in (23),(3)
update the components ofBn+3/2 according to (15),(4)
update the components ofHn+3/2 according to (26).


## 4. Simulation of the Elliptic-Cylindrical Cloak

We calculate electromagnetic propagation for the elliptic-cylindrical cloak using the FDTD formulation shown in Section3. Figure 2 shows the simulation setup: the computational domain is terminated with a perfectly matched layer in the x-direction, and a periodic boundary condition in the y-direction [29]; the inside of the cloak is covered with a perfect electric conductor (PEC); a plane wave source of wavelength λ0= 750 nm (400 THz) is in the TM polarization; the grid spacing is h= 10 nm (λ0/h= 75); and the time step is given by the CFL limit, Δt=h/(c02). Simulation parameters are listed in Table 1.Table 1
Simulation parameters for the elliptic-cylindrical cloak.
Incident wavelengthλ0750 nmInner short semi-axisa500 nmOuter short semi-axisb1000 nmAxis ratiok2Grid spacingh10 nmComputational domain6 × 6μm2Figure 2
Simulation setup of the elliptic-cylindrical cloak for incident TM polarization. Inside of the cloak is covered with a perfect electric conductor (PEC). The gray elliptic region represents the cloak (h: grid spacing).Figure3 shows the FDTD results for the elliptic-cylindrical cloak at the steady state (50 wave periods). Figure 3(a) shows calculated Hz field distributions using h= 10 nm. The wave propagates without significant disturbance around the cloak, and the calculation is stable. The small ripples on phase planes are purely numerical errors and can be made to vanish by reducing the grid spacing. This can be confirmed by calculating the radar cross section (RCS) [29, 51]. In two dimensions, the RCS is defined by
(27)σ(ϕ)=lim⁡r→∞2πr|Es(ϕ)|2|E0|2,
where ϕ is the scattering angle, |Es(ϕ)|2 is the scattered power in far field, and |E0|2 is the incident power. If there is no significant disturbance by the object, σ approaches zero. Figure 3(b) shows normalized RCSs on dB scale, σ/λ0, scattered by a PEC (without cloak) and cloak using different grid spacings, h= 20, 10, and 5 nm. The PEC or cloak using a coarse grid spacing scatters strongly, but the RCS of the cloak rapidly decreases as the grid spacing is reduced.FDTD results for the elliptic-cylindrical cloak. (a)Hz field distributions using h= 10 nm (black lines are inside and outside shells of the cloak). (b) Radar cross section (RCS) scattered by a PEC (without cloak) and cloak using different grid spacings, h= 20, 10, and 5 nm (radial coordinate is dB).
(a)(b)Finally, we examine the cloaking performance of the elliptic-cylindrical cloak using the Drude model. Simulation parameters are the same as Table1 and the cloak is optimized to a wavelength of 750 nm. In the wavelength band, 600 nm–900 nm, we calculate the total cross section (TCS) defined by(28)σt=∫σ(ϕ)dϕ.
Figure 4(a) shows the calculated TCS spectrum. The TCS rapidly increases with wavelength shifts off the optimal. Figure 4(b) shows the RCS for several wavelengths, A: 730 nm, B: 750 nm, and C: 830 nm (normalized to the RCS for 750 nm). For wavelengths A and C, the scattering is much stronger than the optimal wavelength B.Cloaking performance of the elliptic-cylindrical cloak using the Drude model. (a) Total cross section (TCS) spectrum. (b) Sampled RCSs of incident wavelengths, A: 730 nm, B: 750 nm, and C: 830 nm.
(a)(b)


## 5. Conclusion

We describe a stable FDTD modeling procedure for a cloak with a nondiagonal permittivity tensor. When the eigenvalues of the material parameters are less than one, they must be mapped to a dispersion model in order to maintain numerical stability. We implement our method for an elliptic-cylindrical cloak in the TM mode. Numerical calculations demonstrated stable results and the cloaking performance.



---

*Source: 536209-2012-06-04.xml*