Finite-difference time-domain simulation of spacetime cloak

Jason Cornelius; Jinjie Liu; Moysey Brio

doi:10.1364/OE.22.012087

1. Introduction

In recent years, Transformation Optics (TO) has become one of more interesting topics in science. Utilizing transformation optics and metamaterials, spatial invisibility cloaks [1, 2] have been developed. A spatial cloak functions by manipulating the spatial path of electromagnetic waves within the cloaked region in such a way that waves are bent around the object being cloaked. Transformation optics allows for a description of the medium (ie: material parameters) needed for the construction of such devices. For this particular type of cloak, the implementation requires the use of inhomogeneous media where the inhomogeneity is due to the use of spatially dependent anisotropic permeability and permittivity.

Within the past three years, a new type of cloak has been developed [3, 4], the spacetime cloak. Contrary to the spatial cloak, a spacetime cloak functions by altering the speed of the electromagnetic wave as it propagates through the cloaked region. Initially, the speed of the front of the wave is increased while the speed of the back of the wave is decreased. This allows for a gap to occur in which an event can be carried out undetected. To close the gap, the back of the wave is sped up and the front is slowed down allowing the wave to continue on at its original speed. The design of such a cloak can be achieved by transformation optics, and it results in a magneto-electric medium with time and space varying permittivity, permeability, and coupling coefficient.

Due to the interesting nature of these devices, simulation is desirable. To this end, we turn to the Finite-Difference Time-Domain (FDTD) method. The FDTD method [5–7] is a very successful method for simulating electromagnetic wave propagation. It has been applied to the study of various types of materials [7, 8], including dielectrics, linear dispersive materials, nonlinear Raman and Kerr materials, nonlinear dispersive materials [9], bi-isotropic media [10, 11], etc. In this paper, we are interested in space and time dependent magneto-electric materials and their direct application to the spacetime cloak. In [4], the FDTD method has been applied to simulate the spacetime cloak but the authors point out that there are some difficulties near the closing process of the spacetime cloak. Our numerical simulations demonstrate that there are instabilities for conventional FDTD simulation due to the temporal extrapolation of the time dependent magneto-electric constitutive equations.

To stably solve Maxwell’s equations in a time dependent magneto-electric medium, we propose a modified FDTD method based on the use of time overlapped grids. Our method uses two Yee grids that are offset in time by a half time step so that collocated fields are provided to the time dependent magneto-electric constitutive equations. The resulting numerical update equations do not require time extrapolation, therefore the instability due to time extrapolation is avoided. Our method is applied to simulate the spacetime cloak.

2. Transformation optics and spacetime cloak

In this section, we briefly review the derivation of the spacetime cloak using transformation optics. Consider the following covariant form of the Maxwell’s equations (Ampere’s Law and ` the Gauss’ Law of electric fields)

\tilde{\nabla} \tilde{M} = \tilde{\nabla} (\begin{matrix} 0 & D_{x} & D_{y} & D_{z} \\ - D_{x} & 0 & H_{z} & - H_{y} \\ - D_{y} & - H_{z} & 0 & H_{x} \\ - D_{z} & H_{y} & - H_{x} & 0 \end{matrix}) = 0,

where

\tilde{\nabla} = (\frac{\partial}{\partial t}, \frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z})

. Applying a spacetime transformation from (t, x) to (τ, ξ) and keeping y and z coordinates unchanged, we have the following covariant form of the Maxwell’s equation in the new coordinates (τ, ξ, η, ζ):

{\tilde{\nabla}}^{'} {\tilde{M}}^{'} = {\tilde{\nabla}}^{'} (\begin{matrix} 0 & D_{ξ} & D_{η} & D_{ζ} \\ - D_{ξ} & 0 & H_{ζ} & - H_{η} \\ - D_{η} & - H_{ζ} & 0 & H_{ξ} \\ - D_{ζ} & H_{η} & - H_{ξ} & 0 \end{matrix}) = 0,

where

{\tilde{\nabla}}^{'} = (\frac{\partial}{\partial τ}, \frac{\partial}{\partial ξ}, \frac{\partial}{\partial η}, \frac{\partial}{\partial ζ})

, M̃ = |Λ̃|Λ̃⁻¹M̃′Λ̃^−T, and Λ̃ is the Jacobian matrix (I₂ represents the 2 × 2 identity matrix):

\tilde{Λ} = \frac{\partial (τ, ξ, η, ζ)}{\partial (t, x, y, z)} = (\begin{matrix} τ_{t} & τ_{x} & 0 & 0 \\ ξ_{t} & ξ_{x} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}) = (\begin{matrix} Λ & 0 \\ 0 & I_{2} \end{matrix}) .

Similar matrix equations are obtained for Faraday’s Law together with the Gauss’ Law of magnetic fields.

We consider a spacetime cloak design similar to the spacetime cloaks proposed in [3, 4]. Figure 1 shows a diamond shaped spacetime cloak in the (t, x) domain and the transformed (τ, ξ) free space. The original and the transformed regions are the same parallelograms centered at (t₀, x₀) so we have ξ₀ = x₀, τ₀ = t₀, τ_p = t_p, and τ_q = t_q. In the parallelogram we apply the following coordinate transformation to create such a spacetime cloak:

τ = t,

ξ = \frac{σ (x - x_{S})}{σ - δ (t - t_{S}) / (t_{0} - t_{S})} + x_{S},

where the values of x_S and t_S depend on the region where (t, x) lies in:

t_{S} = {\begin{array}{l} t_{p}, & if (t, x) lies in region I or II, \\ t_{q}, & if (t, x) lies in region III or IV, \end{array}

and

x_{S} = {\begin{array}{l} x_{0} + (t - t_{0}) c + σ, & if (t, x) lies in region I or IV, \\ x_{0} + (t - t_{0}) c - σ, & if (t, x) lies in region II or III . \end{array}

Fig. 1 (a) A spacetime cloak in the (t, x) domain. (b) Free space after the transformation to the (τ, ξ) domain.

Download Full Size | PDF

Since the transformed (τ, ξ) domain is free space, we have the constitutive equations D = εE and B = μH. The corresponding constitutive equations in (t, x) domain are magneto-electric:

E_{z} = α D_{z} + β B_{y},

H_{y} = β D_{z} + γ B_{y},

where α, β, and γ are space and time dependent and

{\begin{array}{l} α & = & 1 / a_{12}, \\ β & = & - a_{22} / a_{12}, \\ γ & = & (a_{12} a_{21} - a_{11} a_{22}) / a_{12}, \end{array}

and

(\begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix}) = Λ^{- 1} (\begin{matrix} 0 & ε \\ 1 / μ & 0 \end{matrix}) Λ .

3. The conventional FDTD method for time dependent magneto-electric medium

Consider the following 1+1 D Maxwell’s equations

\frac{\partial D_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x},

\frac{\partial B_{y}}{\partial t} = \frac{\partial E_{z}}{\partial x},

where the constitutive equations are given in Eqs. (8) and (9). For simplicity, we let E = E_z, D = D_z, H = H_y, and B = B_y. We note that the primary issue with simulating magneto-electric media arises from the more complicated constitutive relations given in Eqs. (8) and (9). Using the conventional FDTD method, the discretized Maxwell’s equations (the D first scheme) are

D_{i}^{n + \frac{1}{2}} = D_{i}^{n - \frac{1}{2}} + \frac{Δ t}{Δ x} (H_{i + \frac{1}{2}}^{n} - H_{i - \frac{1}{2}}^{n}),

B_{i + \frac{1}{2}}^{n + 1} = B_{i + \frac{1}{2}}^{n} + \frac{Δ t}{Δ x} (E_{i + 1}^{n + \frac{1}{2}} - E_{i}^{n + \frac{1}{2}}),

and the constitutive relations (assuming that the coupled constitutive relations respond locally in space and time) become

E_{i}^{n + \frac{1}{2}} = α_{i}^{n + \frac{1}{2}} D_{i}^{n + \frac{1}{2}} + β_{i}^{n + \frac{1}{2}} B_{i}^{n + \frac{1}{2}},

H_{i + \frac{1}{2}}^{n + 1} = β_{i + \frac{1}{2}}^{n + 1} D_{i + \frac{1}{2}}^{n + 1} + γ_{i + \frac{1}{2}}^{n + 1} B_{i + \frac{1}{2}}^{n + 1} .

Parameters α, β, and γ are computed analytically from Eq. (10) and are thus available at any point (t, x) in the computational domain. However, the quantities $B_{i}^{n + \frac{1}{2}}$ in Eq. (16) and $D_{i + \frac{1}{2}}^{n + 1}$ in Eq. (17) are not computed in the discretized Maxwell’s equations, Eqs. (14) and (15). As a result, the conventional Yee approach cannot be directly used in the simulation of magneto-electric media.

One method for resolving these issues is to use time extrapolation and space interpolation. To compute $B_{i}^{n + \frac{1}{2}}$ in Eq. (16), we have,

B_{i}^{n + \frac{1}{2}} \approx \frac{3 B_{i}^{n} - B_{i}^{n - 1}}{2}

\approx \frac{\frac{3}{2} (B_{i + \frac{1}{2}}^{n} + B_{i - \frac{1}{2}}^{n}) - \frac{1}{2} (B_{i + \frac{1}{2}}^{n - 1} + B_{i - \frac{1}{2}}^{n - 1})}{2}

\approx \frac{3 B_{i + \frac{1}{2}}^{n} + 3 B_{i - \frac{1}{2}}^{n} - B_{i + \frac{1}{2}}^{n - 1} - B_{i - \frac{1}{2}}^{n - 1}}{4},

which is valid after the first timestep. For the initial timestep, we use

B_{i}^{n}

as an approximation for

B_{i}^{n + \frac{1}{2}}

. Similarly, we compute

D_{i + \frac{i}{2}}^{n + 1}

in Eq. (17) as follows

D_{i + \frac{1}{2}}^{n + 1} \approx \frac{3 D_{i + 1}^{n + \frac{1}{2}} + 3 D_{i}^{n + \frac{1}{2}} - D_{i + 1}^{n - \frac{1}{2}} - D_{i}^{n - \frac{1}{2}}}{4} .

Numerical simulation of the spacetime cloak shows that the conventional FDTD method presented in this section is unstable. This is primarily due to the use of time extrapolation. To avoid such extrapolation, we propose an overlapping Yee algorithm for time dependent magneto-electric media as described in the next section.

4. The overlapping Yee algorithm for time dependent magneto-electric media

Rather than a single D first scheme or B first scheme, we use a combination of both. While computationally more expensive, this produces D and B values at every half time step and spatial step within our computational domain. As a result, we are able to perform the required magneto-electric Maxwell constitutive relation updates without extrapolation.

Figure 2 illustrates the fundamental differences between utilizing a conventional FDTD approach and the overlapping grid approach. Figure 2(a) depicts that, due to the coupled constitutive relations, at any given step, E_z has a dependence on the quantity $\bar{B_{y}}$ . The dotted arrows leading to this value indicate that $\bar{B_{y}}$ must be extrapolated in time and interpolated in space from values known on the conventional FDTD grid.

Fig. 2 Computation of E_z using (a) the time extrapolation under a conventional FDTD method and (b) the overlapping Yee algorithm without extrapolation.

Download Full Size | PDF

On the other hand, Fig. 2(b) illustrates that by utilizing a dual overlapping grid approach, the dependence on time extrapolation and spatial interpolation can be eliminated. A new superscript notation is used to indicate the usage of two overlapping grids and to specify which grid each value of B_y and E_z comes from. All values needed to compute E_z at a particular instance in time exists on one of the two grids.

In presenting an algorithm for the adjusted scheme we utilize the convention that n represents the time coordinate of the grid point and i represents the spatial coordinate. Indices of $n \pm \frac{1}{2}$ and $i \pm \frac{1}{2}$ are use to represent locations that are a half time step or spatial step from the gridpoint.

With this established, we present an algorithm for an overlapping Yee FDTD method.

The algorithm in 1D case can be summarized by the following steps:

Update $D_{i}^{n + \frac{1}{2}}$ on the first Yee grid and $B_{i}^{n + \frac{1}{2}}$ on the second Yee grid: $D_{i}^{n + \frac{1}{2}} = D_{i}^{n - \frac{1}{2}} + \frac{Δ t}{Δ x} [H_{i + \frac{1}{2}}^{n} - H_{i - \frac{1}{2}}^{n}],$ $B_{i}^{n + \frac{1}{2}} = B_{i}^{n - \frac{1}{2}} + \frac{Δ t}{Δ x} [E_{i + \frac{1}{2}}^{n} - E_{i - \frac{1}{2}}^{n}] .$
Update $E_{i}^{n + \frac{1}{2}}$ and $H_{i}^{n + \frac{1}{2}}$ : $E_{i}^{n + \frac{1}{2}} = α_{i}^{n + \frac{1}{2}} D_{i}^{n + \frac{1}{2}} + β_{i}^{n + \frac{1}{2}} B_{i}^{n + \frac{1}{2}},$ $H_{i}^{n + \frac{1}{2}} = β_{i}^{n + \frac{1}{2}} D_{i}^{n + \frac{1}{2}} + γ_{i}^{n + \frac{1}{2}} B_{i}^{n + \frac{1}{2}} .$
Update $B_{i + \frac{1}{2}}^{n + 1}$ on the first Yee grid and $D_{i + \frac{1}{2}}^{n + 1}$ on the second Yee grid: $B_{i + \frac{1}{2}}^{n + 1} = B_{i + \frac{1}{2}}^{n} + \frac{Δ t}{Δ x} [E_{i + 1}^{n + \frac{1}{2}} - E_{i}^{n + \frac{1}{2}}],$ $D_{i + \frac{1}{2}}^{n + 1} = D_{i + \frac{1}{2}}^{n} + \frac{Δ t}{Δ x} [H_{i + 1}^{n + \frac{1}{2}} - H_{i}^{n + \frac{1}{2}}] .$
Update $E_{i + \frac{1}{2}}^{n + 1}$ and $H_{i + \frac{1}{2}}^{n + 1}$ : $E_{i + \frac{1}{2}}^{n + 1} = α_{i + \frac{1}{2}}^{n + 1} D_{i + \frac{1}{2}}^{n + 1} + β_{i + \frac{1}{2}}^{n + 1} B_{i + \frac{1}{2}}^{n + 1},$ $H_{i + \frac{1}{2}}^{n + 1} = β_{i + \frac{1}{2}}^{n + 1} D_{i + \frac{1}{2}}^{n + 1} + γ_{i + \frac{1}{2}}^{n + 1} B_{i + \frac{1}{2}}^{n + 1} .$

We note that the extension to 2-D and 3-D should be feasible by utilizing a time collocated grid (where B and E are computed at the same instance in time, but not necessarily in space). Doing so eliminates time extrapolation in 2-D and 3-D but still uses interpolation in space.

5. Spacetime cloak simulation

The simulation of magneto-electric media has direct applications to the spacetime cloak presented in [3]. We create a similar diamond shaped spacetime cloak using the transformation depicted in Fig. 1. In simulating the above cloak, we utilize a 10 μm by 15 fs long computational domain with a plane wave (with λ = 600 nm) traveling in the positive x direction. Our spacetime cloak is 1200 nm wide and lasts 6 fs. With regards to Fig. 1(a) we have x₀ = 1.5 μm, t₀ = 9 fs, t_q = 6 fs, t_p = 12 fs, σ = 600 nm and δ = 180 nm. The computational domain is terminated by perfectly matched layer (PML) absorbing boundary conditions [7]. Due to the fast and slow phase velocities in the spacetime cloak [3], Δt is chosen by Δt ≤ Δx/v_max where v_max is the largest wave speed computed in the cloaking region. We perform the simulation 3 times with the number of spatial cells, N, equal to 1600, 2000, and 2400 and the same number of time steps

The numerical simulation results are shown in Fig. 3. The figures shown are contour plots of the intensity of the electric field generated by the interaction between the incident source and the spacetime cloak. We compare the conventional FDTD method with time extrapolation and the proposed overlapping Yee FDTD method for various spatial resolutions. The result shows the bending of the electric field (and consequently the electromagnetic wave) in space and time around the cloaked event. The shown simulation results verify that an extrapolation approach is unstable (near the closing process of the cloaked event) and the proposed overlapping FDTD approach provides stable results.

Fig. 3 Contour plots of the electric field intensity for spacetime cloak simulations. Left column: the overlapping Yee FDTD results with grid sizes (a) N = 1600, (c) N = 2000, and (e) N = 2400, respectively. Right column: the corresponding results using a conventional time extrapolation based FDTD method.

Download Full Size | PDF

We now examine the behavior of the electric field as the wave propagates through the diamond shaped space-time region under the same simulation parameters. As the wave enters the cloak it begins to split as the front portion of the wave speeds up and the back portion of the wave slows down. This effect is made possible due to the space-time dependence of the media within the cloaking region. As this separation occurs, it is mirrored in the electric field producing a gap in the electric field. The early development of this gap is shown in Fig. 4(a) near x = 1000 nm. As the wave progresses through the media, the cloak widens in the x direction allowing for the size of this gap to increase as time progresses. The gap in the electric field, as shown in Fig. 4(b), is at a maximum at the center of the cloaking region. As the wave propagates past this point, the cloak begins to close (shrink in the x direction). Due to the closing of the cloak, the front of the wave begins to slow down while the back speeds up. This action allows for the wave to remerge just as it exits the cloak. The field behavior moments before the wave exits the cloak is shown in Fig. 4(c). At this point a small gap near x = 2200 nm still exists. Once the wave has propagated through the cloak region, it is fully rejoined and the propagation behavior and resulting field behavior returns to that of propagation in an isotropic media, as shown in Fig. 4(d).

Fig. 4 Evolution of the Electric field as the wave propagates through the spacetime cloak.

Download Full Size | PDF

6. Stability analysis

To conduct a stability analysis, we compare the distribution of eigenvalues of the conventional FDTD method to that of our overlapping grid approach with respect to the unit circle. For both tests we utilize similar simulation parameters as in the previous section, but with a smaller mesh size of N = 200. We compute the eigenvalues of the numerical update equations at a time instant when t ∈ (t_q, t_p) and plot them in the complex plane. The results when t = 8.34 fs are shown in Fig. 5.

Fig. 5 Distribution of the eigenvalues in the complex plane for (a) the overlapping Yee FDTD method and (b) the conventional FDTD method.

Download Full Size | PDF

As indicated in the Fig. 5(a) all the eigenvalues of the overlapping approach lie within or on the unit circle with the ones lying on unit circle being simple eigenvalues. On the other hand, Fig. 5(b) shows that some of eigenvalues (near the upper and lower right corners) of the conventional FDTD approach lie outside the unit circle. In particular the modulus of the largest eigenvalue is 1 for the overlapping method and 1.046 for the conventional FDTD method. This indicates that the overlapping Yee FDTD method is stable while the conventional FDTD approach (with extrapolation) would experience instability.

7. Conclusion

In conclusion, we have proposed a stable FDTD method for the simulation of space and time dependent magneto-electric medium based on the use of two sets of overlapping Yee grids that are offset in time by a half time step. Additionally, we have shown the direct application of this method to the simulation of the spacetime cloak. The proposed method is useful in exploration of other new physical possibilities offered by metamaterials and transformation optics.

Acknowledgments

We thank Dr. Yong Zeng for invaluable discussions. This work was supported in part by the US Air Force Office of Scientific Research Grant FA9550-10-1-0127, US ARO Grant W911NF-11-2-0046, and NSF Grant HRD-1242067. M. Brio was also supported by US AFOSR MURI Grant FA9550-10-1-0561.

References and links

1. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science 312, 1780–1782 (2006). [CrossRef] [PubMed]

2. U. Leonhardt, “Optical conformal mapping,” Science 312, 1777–1779 (2006). [CrossRef] [PubMed]

3. M. W. McCall, A. Favaro, P. Kinsler, and A. Boardman, “A spacetime cloak, or a history editor,” J. Opt. 13, 024003 (2011). [CrossRef]

4. P. Kinsler and M. W. McCall, “Cloaks, editors, and bubbles: applications of spacetime transformation theory,” Ann. Phys. 526, 51–62 (2014). [CrossRef]

5. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). [CrossRef]

6. A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans. Microwave Theory Tech. 23, 623–630 (1975). [CrossRef]

7. A. Taflove and S. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, Norwood, MA, 2005), 3rd edition.

8. F. L. Teixeira, “Time-domain finite-difference and finite-element methods for Maxwell equations in complex media,” IEEE Trans. Antennas Propag. 56, 2150–2166 (2008). [CrossRef]

9. J. Liu, M. Brio, Y. Zeng, A. Zakharian, W. Hoyer, S. W. Koch, and J. V. Moloney, “Generalization of the FDTD algorithm for simulations of hydrodynamic nonlinear Drude model,” J. Comput. Phys. 229, 5921–5932 (2010). [CrossRef]

10. A. Akyurtlu and D. H. Werner, “BI-FDTD: A novel finite-difference time-domain formulation for modeling wave propagation in bi-isotropic media,” IEEE Trans. Antennas Propag. 52, 416–425 (2004). [CrossRef]

11. A. Semichaevsky, A. Akyurtlu, D. Kern, D. H. Werner, and M. G. Bray, “Novel BI-FDTD approach for the analysis of chiral cylinders and spheres,” IEEE Trans. Antennas Propag. 54, 925–932 (2006). [CrossRef]

Finite-difference time-domain simulation of spacetime cloak

Abstract

1. Introduction

2. Transformation optics and spacetime cloak

3. The conventional FDTD method for time dependent magneto-electric medium

4. The overlapping Yee algorithm for time dependent magneto-electric media

5. Spacetime cloak simulation

6. Stability analysis

7. Conclusion

Acknowledgments

References and links

Cited By

Figures (5)

Equations (29)

Optics Express