FDTD Modeling of a Cloak with a Nondiagonal Permittivity Tensor

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.

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) where 𝜀^𝑖𝑗 is the relative permittivity, 𝜇^𝑖𝑗 is the relative permeability, 𝑔^𝑖𝑗 is the metric tensor, and 𝑔=det𝑔^𝑖𝑗. Because 𝜀^𝑖𝑗,𝜇^𝑖𝑗 are constructed from the symmetric metric tensor 𝑔^𝑖𝑗, they are symmetric. Consequently, 𝜀^𝑖𝑗,𝜇^𝑖𝑗 have real eigenvalues with orthogonal eigenvectors and are thus diagonalizable. The eigenvalues, 𝜆, of 𝜀^𝑖𝑗,𝜇^𝑖𝑗 for an eigenvector 𝐕 are defined by𝜀^𝑖𝑗𝐕=𝜇^𝑖𝑗𝐕=𝜆𝐕.(2) The phase velocity of light in a material is given by 𝑐=𝑐₀/𝜆 (𝑐₀= vacuum light speed), and the Courant-Friedrichs-Lewy (CFL) stability limit becomesΔ𝑡≤𝜆ℎ𝑐₀√𝑑,(3) where Δ𝑡 is the time step, ℎ is the grid spacing, and 𝑑= 1, 2, and 3 dimensions. Since the FDTD stability depends on the eigenvalues of 𝜀^𝑖𝑗 and 𝜇^𝑖𝑗, to analyze nondiagonal cases, we must first find the eigenvalues and diagonalize 𝜀^𝑖𝑗 and 𝜇^𝑖𝑗. After the diagonalization, the FDTD method for diagonal cases [31–38] can be applied. For diagonal 𝜀^𝑖𝑗 and 𝜇^𝑖𝑗, 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 𝜀^𝑖𝑗 and 𝜇^𝑖𝑗 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 𝜀^𝑖𝑗 and 𝜇^𝑖𝑗 in orthonormal elliptic-cylindrical coordinates [45, 46], and in the other 𝜀^𝑖𝑗 and 𝜇^𝑖𝑗 are nondiagonal in Cartesian coordinates [47, 48]. We derive a FDTD formulation for the latter in the transverse magnetic (TM) polarization.

3.1. Diagonalization

Figure 1 shows an elliptic-cylindrical cloak in Figure 1(a) Cartesian coordinates and Figure 1(b) transformed coordinates. The inner axis 𝑎, the outer axis 𝑏, and the perpendicular axes 𝑘𝑎 and 𝑘𝑏 are depicted. The elliptic-cylindrical cloak is horizontal when 𝑘>1, and vertical when 𝑘<1. In the cloak region, 𝑘𝑎≤√𝑥²+𝑘²𝑦²≤𝑘𝑏, the material parameters are expressed by𝜀^𝑖𝑗=𝜇^𝑖𝑗=⎡⎢⎢⎢⎢⎣𝜀^𝑥𝑥𝜀_𝑥𝑦𝜀_𝑥𝑦00𝜀0_𝑦𝑦𝜀0_𝑧𝑧⎤⎥⎥⎥⎥⎦,(4) where 𝜀_𝑥𝑥=𝑟𝑟−𝑘𝑎+𝑘²𝑎²𝑅²−2𝑘𝑎𝑟³(𝑟−𝑘𝑎)𝑟⁵𝜀_𝑥𝑦𝑥²,=𝑘²𝑎²𝑅²−𝑘𝑎1+𝑘²𝑟³(𝑟−𝑘𝑎)𝑟⁵𝜀𝑥𝑦,_𝑦𝑦=𝑟𝑟−𝑘𝑎+𝑘²𝑎²𝑅²−2𝑘³𝑎𝑟³(𝑟−𝑘𝑎)𝑟⁵𝜀𝑦²,(5)_𝑧𝑧=𝑏𝑏−𝑎²𝑟−𝑘𝑎𝑟,(6) where 𝑟=√𝑥²+𝑘²𝑦² and 𝑅=√𝑥²+𝑘⁴𝑦². From (4) to (6) we can obtain three eigenvalues 𝜆₁=𝛼−1𝛼+1,𝜆²=1𝜆₁,𝜆₃=𝜀_𝑧𝑧,(7) where 𝛼=1+4𝑟⁵(𝑟−𝑘𝑎)𝑘²𝑎²𝑅²𝑥²+𝑦².(8) Since 𝜀^𝑖𝑗 is symmetric, it is diagonalized by the eigenvalue matrix Λ and its orthogonal matrix 𝑃 as follows: 𝜀^𝑖𝑗=𝑃Λ𝑃^𝑇,(9) where Λ=⎡⎢⎢⎢⎢⎣𝜆¹0𝜆0₂000𝜆0₃𝑃=⎡⎢⎢⎢⎢⎢⎢⎣𝜀_𝑥𝑦⎤⎥⎥⎥⎥⎦,(10)𝛽𝜆²−𝜀_𝑦𝑦𝜀_𝑦𝑦𝛽0−𝜆₂𝛽0𝜀𝛽0^𝑥𝑦10⎤⎥⎥⎥⎥⎥⎥⎦,(11) where 𝛽=𝜀^2𝑥𝑦+(𝜆₂−𝜀_𝑦𝑦)².

(a)

(b)

Figure 1 

Elliptic-cylindrical cloak. (a) Cartesian coordinates: 𝑎, 𝑏 are inner and outer axes; 𝑘𝑎, 𝑘𝑏 are the perpendicular axes. (b) Transformed coordinates.

3.2. Mapping Eigenvalues to a Dispersion Model

From (7) and (8), 𝜆₁ and 𝜆₃ have values less than one in the cloak region (𝑘𝑎≤𝑟≤𝑘𝑏). Thus, 𝜆₁,𝜆₃ must be replaced by dispersive quantities by using (for example) the Drude model 𝜆_𝑖=𝜀_∞𝑖−𝜔^2𝑝𝑖𝜔²−𝑗𝜔𝛾_𝑖,(𝑖=1,3),(12) where 𝜔 is the angular frequency, 𝜀_∞𝑖 is the infinite-frequency permittivity, 𝜔_𝑝𝑖 is the plasma frequency, and 𝛾_𝑖 is the collision frequency. For simplicity, we consider the lossless case, 𝛾_𝑖=0. Then the plasma frequencies are given by 𝜔_𝑝𝑖=𝜔√𝜀_∞𝑖−𝜆_𝑖, where 𝜀_∞𝑖=max(1,𝜆_𝑖).

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) where 𝐃 is the electric flux density, 𝐇 is the magnetic field, 𝐁 is the magnetic flux density, and 𝐄 is the electric field. In the TM polarization, electromagnetic fields reduce to three nonzero components 𝐸_𝑥, 𝐸_𝑦, and 𝐻_𝑧 (𝐷_𝑥, 𝐷_𝑦, and 𝐵_𝑧). The 𝐃- and 𝐁-update equations are obtained using Yee algorithm [26–29] as follows: 𝐷^𝑥𝑛+1=𝐷^𝑛𝑥+Δ𝑡ℎ𝑑^𝑦𝐷_𝑦𝑛+1𝐻^𝑧𝑛+1/2,=𝐷^𝑛𝑦−Δ𝑡ℎ𝑑^𝑥𝐵𝐻^𝑧𝑛+1/2,(14)_𝑧𝑛+3/2=𝐵^𝑧𝑛+1/2−Δ𝑡ℎ𝑑^𝑥𝐸^𝑦𝑛+1−𝑑_𝑦𝐸^𝑥𝑛+1,(15) where we simply write 𝐷_𝑥(𝑡=𝑛Δ𝑡)→𝐷^𝑛𝑥 (𝑛= integer) and 𝑑_𝑥, 𝑑_𝑦 are the spatial difference operators defined by 𝑑_𝑥𝑓(𝑥,𝑦)=𝑓𝑥+ℎ𝑑2,𝑦−𝑓𝑥−ℎ2,𝑦_𝑦,𝑓(𝑥,𝑦)=𝑓𝑥,𝑦+ℎ2−𝑓𝑥,𝑦−ℎ2.(16)

To find the 𝐄-update equations, we consider the relation 𝐃=𝜀₀𝜀^𝑖𝑗𝐄,(17) where 𝜀₀ is the vacuum permittivity. From (9), we obtain 𝜀₀𝐄=𝜀^𝑖𝑗⁻¹𝐃=𝑃Λ⁻¹𝑃^𝑇𝐃.(18) Substituting (10) in (18) and multiplying 𝜆₁𝜆₂ by both sides, we obtain 𝜀₀𝜆₁𝜆₂𝐸_𝑥=𝜆₁𝑡²²+𝜆₂𝑡²¹𝐷_𝑥+𝑡₁𝑡₂𝜆₁−𝜆₂𝐷_𝑦,(19)𝜀₀𝜆₁𝜆₂𝐸_𝑦=𝜆₁𝑡²¹+𝜆₂𝑡²²𝐷_𝑦+𝑡₁𝑡₂𝜆₁−𝜆₂𝐷_𝑥,(20) where 𝑡₁=𝜀_𝑥𝑦/𝛽 and 𝑡₂=(𝜆₂−𝜀_𝑦𝑦)/𝛽. Substituting the Drude model for 𝜆₁ as shown in (12) and using the inverse Fourier transformation rule, −𝜔²→𝜕²/𝜕𝑡², (19) becomes 𝜀₀𝜆₂𝜀^∞1𝜕²𝜕𝑡²+𝜔^2𝑝1𝐸^𝑥=𝜀_∞1𝑡²²+𝜆₂𝑡²¹𝜕²𝜕𝑡²+𝜔^2𝑝1𝑡²²𝐷^𝑥+𝑡₁𝑡₂𝜀^∞1−𝜆₂𝜕²𝜕𝑡²+𝜔^2𝑝1𝐷^𝑦.(21) For the discretization, we use the central difference approximation and the central average operator, 𝜕²𝜕𝑡²𝐸^𝑛𝑥=𝐸^𝑥𝑛+1−2𝐸^𝑛𝑥+𝐸^𝑥𝑛−1Δ𝑡²𝐸_𝑛𝑥,=𝐸^𝑥𝑛+1+2𝐸^𝑛𝑥+𝐸^𝑥𝑛−14.(22) The central average operator improves the stability and accuracy [40, 49, 50]. Similarly, 𝐷_𝑥 and 𝐷_𝑦 are discretized, and we obtain the 𝐸_𝑥-update equation 𝐸^𝑥𝑛+1=−𝐸^𝑥𝑛−1+2𝑎⁻¹𝑎⁺¹𝐸^𝑛𝑥+1𝜀₀𝜆₂𝑎⁺¹×𝑏⁺¹𝐷^𝑥𝑛+1+𝐷^𝑥𝑛−1−2𝑏⁻¹𝐷^𝑛𝑥+𝑐⁺¹𝐷^𝑦𝑛+1+𝐷^𝑦𝑛−1−2𝑐⁻¹𝐷^𝑛𝑦,(23) where 𝐷^𝑦𝑛+1,𝐷^𝑛𝑦,𝐷^𝑦𝑛−1 must be spatially interpolated due to the staggered Yee cell [31], and 𝑎^±𝑖=𝜀^∞𝑖Δ𝑡²±𝜔^2𝑝𝑖𝑏_±𝑖4,=𝜀^∞𝑖𝑡²²+𝜆₂𝑡²¹Δ𝑡²±𝜔^2𝑝𝑖𝑡²²𝑐_±𝑖4,=𝑡₁𝑡₂𝜀_∞𝑖−𝜆₂Δ𝑡²±𝜔^2𝑝𝑖𝑡²²4.(24) Similarly, the 𝐸_𝑦-update equation is obtained by exchanging 𝑡₁↔𝑡₂, 𝐸_𝑥↔𝐸_𝑦, and 𝐷_𝑥↔𝐷_𝑦.

To find the 𝐻_𝑧-update equation, we consider the relation 𝐵_𝑧=𝜇₀𝜀_𝑧𝑧𝐻_𝑧=𝜇₀𝜀^∞3−𝜔^2𝑝3𝜔²𝐻^𝑧.(25) Analogously to the 𝐸_𝑥 field, the 𝐻_𝑧-update equation can be obtained in the form 𝐻^𝑧𝑛+3/2=−𝐻^{𝑧𝑛−1/2}+2𝑎⁻³𝑎⁺³𝐻^𝑧𝑛+1/2+𝐵^𝑧𝑛+3/2−2𝐵^𝑧𝑛+1/2+𝐵^{𝑧𝑛−1/2}𝜇₀Δ𝑡²𝑎⁺³,(26) where 𝜇₀ is the vacuum permeability.

In summary, the electromagnetic fields are iteratively updated in the following sequence:(1)update the components of 𝐃^𝑛+1 according to (14),(2)update the components of 𝐄^𝑛+1 according to the sample given in (23),(3)update the components of 𝐁^𝑛+3/2 according to (15),(4)update the components of 𝐇^𝑛+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 Section 3. Figure 2 shows the simulation setup: the computational domain is terminated with a perfectly matched layer in the 𝑥-direction, and a periodic boundary condition in the 𝑦-direction [29]; the inside of the cloak is covered with a perfect electric conductor (PEC); a plane wave source of wavelength 𝜆₀= 750 nm (400 THz) is in the TM polarization; the grid spacing is ℎ= 10 nm (𝜆₀/ℎ= 75); and the time step is given by the CFL limit, Δ𝑡=ℎ/(𝑐₀√2). Simulation parameters are listed in Table 1.

Table 1 

Simulation parameters for the elliptic-cylindrical cloak.

Figure 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 (ℎ: grid spacing).

Figure 3 shows the FDTD results for the elliptic-cylindrical cloak at the steady state (50 wave periods). Figure 3(a) shows calculated 𝐻_𝑧 field distributions using ℎ= 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 𝜎(𝜙)=lim_𝑟→∞2𝜋𝑟||𝐸^𝑠(𝜙)||²||𝐸₀||²,(27) where 𝜙 is the scattering angle, |𝐸_𝑠(𝜙)|² is the scattered power in far field, and |𝐸₀|² is the incident power. If there is no significant disturbance by the object, 𝜎 approaches zero. Figure 3(b) shows normalized RCSs on dB scale, 𝜎/𝜆₀, scattered by a PEC (without cloak) and cloak using different grid spacings, ℎ= 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.

(a)

(b)

Figure 3 

FDTD results for the elliptic-cylindrical cloak. (a) 𝐻𝑧 field distributions using ℎ= 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, ℎ= 20, 10, and 5 nm (radial coordinate is dB).

Finally, we examine the cloaking performance of the elliptic-cylindrical cloak using the Drude model. Simulation parameters are the same as Table 1 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) 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.

(a)

(b)

Figure 4 

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.

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.

Acknowledgment

The authors deeply appreciate the financial support of Grant-in-Aid for Japan Society for the Promotion of Science (JSPS) Fellows.