International Journal of Antennas and Propagation
Volume 2020 (2020), Article ID 8294769, 16 pages
https://doi.org/10.1155/2020/8294769

Research Article

Reduction of FDTD Stair-Casing Error regarding to Metamaterials by Using High-Order Polynomial Transformation Functions

Ki-Woong Bae,¹ Dong-Wook Seo,² and Noh-Hoon Myung³

¹Maritime Technology Research Institute, Agency for Defense Development (ADD), Jinhae-gu, Changwon-si, Gyeongsangnam-do, Republic of Korea
²Department of Radio Communication Engineering, Korea Maritime and Ocean University (KMOU), 727 Taejong-ro, Yeongdo-gu, Busan 49112, Republic of Korea
³School of Electrical Engineering, KAIST, 291 Daehak-ro, Yuseong-gu, Daejeon 34141, Republic of Korea

Correspondence should be addressed to Dong-Wook Seo; dwseo@kmou.ac.kr

Received 9 November 2019; Accepted 24 December 2019; Published 13 January 2020

Academic Editor: Giuseppe Castaldi

Copyright © 2020 Ki-Woong Bae et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The material parameters of a metamaterial (MTM) are determined by the transformation function used in the optical transformation. Some previously reported MTMs, such as the invisibility cloak, the field rotator, and the field concentrator, were designed by a linear transformation. Their impedance was matched to the background so that no reflection was found; however, the material parameters were mismatched to the background due to the linear transformation function. In the present work, the parameters were matched by using high-order polynomial functions as the transformation function. Since similar materials are filled in boundary cells of the finite-difference time-domain (FDTD) algorithm, the stair-casing error was reduced and the tolerance against boundary abrasion was increased. The frequency response of the proposed method was analyzed. The proposed method is applicable to MTM structures that have complex boundary shapes. In this work, circular and elliptic boundary shapes were considered as examples.

1. Introduction

The finite-difference time-domain (FDTD) algorithm models an analysis domain on a cell basis; therefore, at a boundary of two distinct media, the FDTD method should handle the boundary with either partially filled cells or simple stair-casing approximation. The partially filled cells require special treatments [1, 2], such as subpixel smoothing [3–5]. The stair-casing approximation fills a boundary cell with one of its filling media; the stair-casing approximation is intuitive and fast at the price of accuracy.

The material parameters of a metamaterial (MTM) are determined by the function used in the optical transformation [6]. Some previously reported MTMs, such as the invisibility cloak, the field rotator, and the field concentrator [7–9], were designed by a linear transformation. In [7–9], their impedance, Z = , was matched to the background; hence, no reflection was found. The material parameters, however, were mismatched to the background. Consequently, stair-case approximating MTM boundaries produced errors in the FDTD algorithm, and for this reason, MTM structures that had complex boundary shape could not be accurately analyzed with FDTD.

The linear function has been utilized in most research regarding MTMs [7–12], and the use of other types of function has not been found often. For example, Luo et al. used functions that are differentiable at boundaries to eliminate scattered fields [13]. Cai et al. used a second-order polynomial to develop a nonmagnetic cloak [14]. Few or no studies, however, have considered the proper type of functions to match material parameters.

In the optical transformation, the derivative of transformation functions should be continuous at boundaries to match material parameters, as will be discussed in detail in Section 2. High-order polynomial functions can be a great help in matching these parameters due to their continuous derivative. The proper type of function will be suggested in Section 3.

As a result of matching parameters, the stair-casing error in the FDTD method can be reduced because similar materials fill the partially filled cells. Mismatched MTMs have the problem of high sensitivity to boundary abrasion. Matching parameters reduces this sensitivity; therefore, the tolerance of MTMs can be enhanced against boundary abrasion.

2. Limitation of Linear Functions

The material parameters of MTMs are calculated by the optical transformation given in (1) and (2) [6]. Here, the space is transformed by a function f. The transformation function f should be continuous at boundary with the space before the transformation to ensure the continuity of space:

where

xi and xi′ represent coordinate components before and after the transformation, respectively. The roman indices i and j run from (1) to (3), for the three spatial components; ε and μ denote material parameters before the transformation, and δ denotes Kronecker’s delta.

If f ′ is continuous at a boundary with the derivative of space before the transformation ∂xi/∂xi, that is,  = I, then parameter matching is achieved as  = ε and  = μ. Here, f ′ is the derivative of f with respect to xi and I is the identity matrix.

The transformation functions of the invisibility cloak, the field rotator, and the field concentrator are cited from [7–9] by (5)–(7). Each of these MTMs has two concentric circular boundaries, which are illustrated in Figures 1(a) and 1(b):

where means the inverse function of f. These linear transformation functions are continuous at both boundaries, R₁ and R₂ (R₁ and R₃ for the concentrator).

(a)

(b)

(c)

Figure 1: (a) Coordinates of the invisibility cloak and the field rotator, (b) coordinates of the field concentrator, and (c) the cells that are replaced with the background material to simulate boundary abrasion.

However, not only the transformation function itself but also its derivative should be continuous at both boundaries to match the material parameters; a total of four conditions are required to be satisfied. A linear function f₁(r) = Ar + B has two coefficients, A and B; thus, the four conditions cannot always be met. On the other hand, an Nth-order polynomial function has N + 1 coefficients. These coefficients can be freely adjusted to meet the required conditions.

3. Proposal of High-Order Polynomial Functions

3.1. Third-Order Polynomial Functions

The aforementioned four conditions can be met by using a third-order polynomial function given as (8). The function has four coefficients, A, B, C, and D. Polynomials that have more coefficients than conditions are not recommended because they vary nonmonotonically and thereupon the material distribution becomes complicated:

where coefficients A, B, C, and D are determined by the substitution of the four boundary conditions about f and f ′ into (5) at both boundaries. The conditions are given as follows. The invisibility cloak is in contact with the background at R₂ and expands the origin into R₁. In mathematical expressions,  = 1 at r′ = R₂ and , respectively. Thus, the conditions should be  = 1 and  = 0, where means . Next, since the rotator does not rotate fields in r′ > R₂ and r′ < R₁, the rotation rate should be 0 at both boundaries, which can be written as  =  = 0. Finally, the concentrator is in contact with the background at R₃ and with the core at R₁, which is expressed as  = 1 and  = R₁/R₂, respectively. These conditions are, respectively, summarized as follows:

The problem corresponds to a system of equations that has four unknowns, A, B, C, and D. These unknowns can be determined by an inverse matrix which has (4 × 4) size. For example, the unknown coefficients of the concentrator can be obtained from the following matrix equation:

The unknown coefficients were derived with the radii given as 3 × R₁ = 1.5 × R₂ = R₃. These are summarized in Table 1.

Table 1: Coefficients of third-order polynomial functions.

Cloak
Rotator
Concentrator

A
−278.083
1784.88
−231.736

B
134.789
−834.249
182.928

C
−19.00
115.533
−45.500

D
0.83093
−4.00
3.73918

The graphs of linear and third-order functions are compared in Figure 2. For example, the linear function f₁ in Figure 2(a) has discontinuous derivatives, such as , while f₃ shows continuous derivatives, such as and . Note that f₃ is monotonic in range 0 < r < R₂; thus, .

(a)

(b)

(c)

Figure 2: Transformation functions: (a) the cloak, (b) the rotator, and (c) the concentrator.

The material parameters of MTMs are calculated by (1) and presented in Figure 3. As a result of the continuous derivative, the parameters with the use of f₃ are continuous, although those with the use of f₁ are discontinuous.

(a)

(b)

(c)

Figure 3: Metamaterial parameters: (a) the cloak, (b) the rotator, and (c) the concentrator.

3.2. Fifth-Order Polynomial Functions

The FDTD stair-casing error originates from abrupt material parameter change; hence, the error can be further reduced by suppressing the parameter change near boundaries. In this section, the material parameters are made to have small change near boundaries to reduce the stair-casing error.

To suppress the parameter change near two boundaries, we propose an additional derivative condition for each boundary, that is, two additional conditions are proposed for each MTM. The MTM in the ranges R₁ < r′ < R₁ +  and R2,3 − < r′ < R2,3 has the possibility of being situated in partially filled cells. Thus, it is logical to set conditions at r′ = R₁ +  and R2,3 − , where Δx denotes the FDTD spatial step size. Six (that is, 4 + 2) conditions can be met by using a fifth-order polynomial function:

If  =  and  = , then f₅ = f₃. If  =  and  = , then the change is completely suppressed because  =  and  = . Thus, it is logical to choose value in the range of and . In (11), 0 ≤ P ≤ 1 is used to express and conveniently. Note that the change is more suppressed with increasing :

The function f₅, on the other hand, is nonmonotonic with high values. This can be explained by Rolle’s theorem and inflection points. A detailed explanation will be omitted for brevity. The upper limits of that make f₅ monotonic were calculated by increasing from 0.01 to 0.99. The calculation reveals that the maximum values are 0.70, 0.55, and 0.63, for the cloak, the rotator, and the concentrator, respectively.

Coefficients A through F were derived from a 6 × 6 matrix similar to (12) with the aforementioned values. These are summarized in Table 2. Figures 4 and 5 reveal that the f₅ and its material parameter change near boundaries are smaller than those obtained with f₃. Notice that the parameter change in the middle is larger with the use of f₅ than that with the use of f3, as shown in Figure 5.

Table 2: Coefficients of fifth-order polynomial functions.

Cloak ()
Rotator ()
Concentrator ()

A
7.6312e + 04
−4.4973e + 05
5.7501e + 04

B
−5.9825e + 04
3.4810e + 05
−7.5034e + 04

C
1.8069e + 04
−1.0356e + 05
3.8622e + 04

D
−2.6118e + 03
1.4726e + 04
−9.7941e + 03

E
1.8148e + 02
−1.0059e + 03
1.2247e + 03

F
−4.8800
27.5782
−60.4071

(a)

(b)

(c)

Figure 4: Proposed fifth-order transformation function: (a) the cloak, (b) the rotator, and (c) the concentrator. The function had smaller changes near boundaries compared with f₃.

(a)

(b)

(c)

Figure 5: Material parameters of the fifth-order function: (a) the cloak, (b) the rotator, and (c) the concentrator. The parameters had smaller changes near boundaries and larger changes in the middle compared to those obtained with f₃.

The stair-casing error in the FDTD algorithm originates from the abrupt material parameter change. Hence, exactly how the stair-casing error was reduced can be shown by visualizing how similar the material parameters are in neighboring cells. Figure 6 depicts the εr component of the field rotator around a stair-cased R₂ boundary. The size of grid is given in Table 3. The εr value changes abruptly at the boundary with the linear function in Figure 6(a), while the change is much smoother with the third-order function in Figure 6(b). The fifth-order function shows very small parameter change across the boundary in Figure 6(c). This tendency was observed in of the three MTMs in consideration. The smoother parameter change indicates that it will produce smaller FDTD stair-casing errors.

(a)

(b)

(c)

Figure 6: The εr component of the field rotator around the R₂ boundary: (a) the linear function, (b) the third-order function, and (c) the fifth-order function. The circular boundary (red dashed line) was approximated to a stair-cased boundary (solid blue line).

Table 3: FDTD simulation conditions.

R₁ (mm)
R₂ (mm)
R₃ (mm)
f₀ (GHz)
∆x (=∆y)

Cloak
106
212
—
2
λ₀/30

Rotator
106
212
—
2
λ₀/30

Concentrator
106
212
318
2
λ₀/30

The proposed high-order function method is applicable to elliptic coordinate systems, as shown in Figure 7(a) [15]. For example, consider a cloak in an elliptic system with a focal length . The values 0.1 and 0.2 m were chosen for the lengths of major axes of inner and outer ellipses, respectively. The linear function (2) in [15] was substituted with the following third-order function similar to (8):

(a)

(b)

(c)

Figure 7: Application to an elliptic coordinate system: (a) the elliptic coordinate, (b) the material parameters obtained with (12), and (c) the μz component obtained with (12) in a stair-cased boundary.

The unknown coefficients were derived from the continuity conditions of f and f ′ as follows: A = –8.47518, B = 44.70059, C = –74.56999, and D = 40.05591.

The material parameters are provided in Figure 7(b). They were finite and changed more smoothly in contrast to those in the circular cloak in Figure 3(a), and most importantly, they were matched to the background. Figure 7(c) illustrates the μz component around a stair-cased outer boundary. The material parameter change was continuous in contrast to that obtained with the linear function.

The linear function produced 0 < μz < 1.88, while the third-order function produced 1.0 < μz < 8.7. The elliptic cloak becomes easier to implement in FDTD with this nondispersive μz component.

4. Numerical Results

The validity of the proposed function was confirmed with FDTD simulations. The 2D FDTD domain is provided in Figure 8 with simulation conditions summarized in Table 3. The temporal step size Δt was adjusted to meet the Courant–Friedrich–Lewy stability condition. The analysis domain was discretized with the second-order FDTD method and truncated with the second-order Mur absorbing boundary condition. The MTMs were stair-casing approximated. The dislocated FDTD fields in the middle were stably interpolated as claimed in [16] to resolve the FDTD instability issue.

Figure 8: Two-dimensional FDTD analysis domain. A source was incident from the left.

The FDTD snapshots with the use of linear and third-order functions are presented in Figures 9 and 10, respectively. The snapshots were taken at the same time step. The error in the form of static charges was previously reported in [17], and the same form of error was found as seen in Figure 9. On the contrary, no stair-casing error was observed in Figure 10. The fields actually changed more smoothly over boundaries in Figure 10 compared with Figure 9. These results are attributable to the matched parameters. The field changed even more smoothly with the use of the fifth-order function, which is attributable to the smaller material parameter change near boundaries.

(a)

(b)

(c)

Figure 9: FDTD snapshots with the use of linear functions: (a) the cloak, (b) the rotator, and (c) the concentrator. The red dotted circles indicate the location of stair-casing error.

(a)

(b)

(c)

Figure 10: FDTD snapshots with the use of third-order functions: (a) the cloak, (b) the rotator, and (c) the concentrator. The field changed smoothly at the boundaries without errors.

The relative root mean square error (RMSE) during the steady state of the simulation was calculated by (16), as shown in Figure 11:

where K is the number of FDTD cells. The error was proportional to since the second-order FDTD algorithm was used [18]. The fifth-order function had the smallest error followed by the third-order and linear functions; the high-order functions reduced the error as intended. On average, the RMSE was reduced about 6.0 dB, 2.5 dB, and 1.4 dB for the cloak, the rotator, and the concentrator, respectively.

(a)

(b)

(c)

Figure 11: Relative root mean square error: (a) the cloak, (b) the rotator, and (c) the concentrator. The fifth-order function had the smallest error followed by the third-order and linear functions.

The cloak shows the largest difference of material parameters at the boundary between the background media followed by the rotator and the concentrator. This explains the steepness of error decrease in Figure 11(a) followed by Figures 11(b) and 11(c).

The reduced error improved calculation accuracy, which resulted in the increase of stable simulation time as provided in Figure 12. The fifth-order function had larger time increase than the third-order function. The maximum stable time was increased about 2.5 times, 3.3 times, and 1.4 times for the cloak, the rotator, and the concentrator, respectively.

(a)

(b)

(c)

Figure 12: Stable simulation time: (a) the cloak, (b) the rotator, and (c) the concentrator. The stable time was increased as the error was reduced.

As seen in Figure 5, the material parameter change was increased in the middle; the FDTD error increases with the decrease of modeling accuracy. Nonetheless, the decrease in RMSE indicates that the amount of error reduction at boundaries was larger than the error increase in the middle. As long as the boundary error is kept smaller than the middle error, these effects are expected to be further enhanced by further suppressing the change near boundaries.

The proposed high-order polynomial functions were applied to lossless cases, but the application to lossy cases is straightforward. In lossy cases, incident energy was dissipated in the MTMs, as the result, the instability might appear in later time steps.

Some boundary cells of the MTMs were replaced with the background media (air) to simulate the boundary abrasion. The replaced cells are marked in Figure 1(c). With the use of the linear function, the scattering pattern of the cloak was similar to that of a sphere, and the fields of the rotator and the concentrator were severely distorted in the core region, as shown in Figure 13. On the contrary, the functionalities remained relatively intact with the use of the third-order function, as shown in Figure 14. These results can be interpreted as indicating the sensitivity is reduced or the tolerance is enhanced against boundary abrasion as a consequence of the matched material parameters.

(a)

(b)

(c)

Figure 13: FDTD snapshots with boundary abrasion with the use of linear functions: (a) the cloak, (b) the rotator, and (c) the concentrator. The MTMs lost their functionality.

(a)

(b)

(c)

Figure 14: FDTD snapshots with boundary abrasion with the use of third-order functions: (a) the cloak, (b) the rotator, and (c) the concentrator. The functionalities of MTMs remained relatively intact.

Frequency responses were analyzed. A broadband Gaussian pulse was launched with a center frequency at 2 GHz and a full width at half maximum (FWHM) bandwidth of 1 GHz. The Hz field was spatially averaged along a line parallel to the y-axis. The averaged field was Fourier transformed and then divided by the input Gaussian spectrum to obtain the transmitted field amplitude; the results are provided in Figure 15. The scattering amplitude of a PEC cylinder with radius R₁ is also included in the graph.

(a)

(b)

(c)

Figure 15: The frequency response of the various functions: (a) the invisibility cloak, (b) the field rotator, and (c) the field concentrator. The insets show the bandwidth characteristics clearly.

The bandwidth of a transformation device was defined as the range of frequencies, where the transmitted amplitude is larger than that of a bare PEC cylinder [19]. According to [19], the bandwidth of an MTM device decreases with the thickness of the device. A higher-order function has a larger material parameter change in the middle (see Figure 5), namely, narrower actual device thickness. This indicates that the linear function had the largest bandwidth followed by the third-order and fifth-order functions. The higher transmitted amplitude in Figure 15(a) indicates that the higher-order functions improved the performance of the invisibility cloak at its design frequency. Despite the actually narrower device thickness, the bandwidth of the rotator and concentrator were merely unaffected.

The stair-casing error in the FDTD method and the sensitivity at boundaries increase with the complexity of boundary shape. The proposed method lowers the material contrast at boundaries, so the definition of boundary shape is decreased. Therefore, the method would be effective for structures that have complex boundary shapes.

In this work, no averaging or interpolation scheme was considered for boundary cells. However, a spatial average at boundaries might improve the FDTD analysis accuracy of the proposed high-order functions [20]. The proposed method can be used to help the design of metasurfaces [21, 22]. The present work focused on the frequency nonlinearity of the MTMs. The field intensity nonlinearity can be studied further in the next paper.

5. Conclusion

Previous studies about MTMs have focused on developing exotic functionalities with the use of linear functions. In the present study, useful properties, such as matched material parameters and impedance were added to MTMs by using high-order polynomials. As a result, the FDTD stair-casing error and the sensitivity at circular and elliptic boundaries were reduced and the stable simulation time was increased. This is the first study about reduction of FDTD error and sensitivity at MTM boundaries. These effects can be enhanced by using higher-order polynomials that have smaller material change near boundaries. These results imply that the proposed method would permit accurate FDTD analysis of complex-shaped MTM structures and would enable the design of such structures with boundary tolerance. The bandwidth characteristic was merely unaffected by the proposed method.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Research Foundation of Korea (NRF) and funded by the Korea government (MSIT) (no. 2018R1C1B6003854).

References

M. Ahmad, N. Hamid, and A. Mario, “Contour-path effective permittivities for the two-dimensional finite-difference time-domain method,” Optics Express, vol. 13, no. 25, pp. 10367–10381, 2005.

J. Nadobny, D. Sullivan, and P. Wust, “A general three-dimensional tensor FDTD-formulation for electrically inhomogeneous lossy media using the Z-transform,” IEEE Transactions on Antennas and Propagation, vol. 56, pp. 1027–1040, 2008.

A. F. Oskooi, C. Kottke, and S. G. Johnson, “Accurate finite-difference time-domain simulation of anisotropic media by subpixel smoothing,” Optics Letters, vol. 34, no. 18, pp. 2778–2780, 2009.

A. Deinega and I. Valuev, “Subpixel smoothing for conductive and dispersive media in the finite-difference time-domain method,” Optics Letters, vol. 32, no. 23, pp. 3429–3431, 2007.

J. Liu, M. Brio, and J. V. Moloney, “Subpixel smoothing finite-difference time-domain method for material interface between dielectric and dispersive media,” Optics Letters, vol. 37, no. 22, pp. 4802–4804, 2012.

D. Schurig, J. B. Pendry, and D. R. Smith, “Calculation of material properties and ray tracing in transformation media,” Optics Express, vol. 14, no. 21, pp. 9794–9804, 2006.

J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science, vol. 312, no. 5781, pp. 1780–1782, 2006.

H. Chen and C. T. Chan, “Transformation media that rotate electromagnetic fields,” Applied Physics Letters, vol. 90, no. 24, p. 241105, 2007.

W. X. Jiang et al., “Design of arbitrarily shaped concentrators based on conformally optical transformation of nonuniform rational B-spline surfaces,” Applied Physics Letters, vol. 92, no. 26, p. 264101, 2008.

Y. Cui, C. Argyropoulos, and Y. Hao, “Full-wave finite-difference time-domain simulation of electromagnetic cloaking structures,” Optics Express, vol. 16, no. 9, pp. 6717–6730, 2008.

C. Argyropoulos, Y. Zhao, and Y. Hao, “A radially-dependent dispersive finite-difference time-domain method for the evaluation of electromagnetic cloaks,” IEEE Transactions on Antennas and Propagation, vol. 57, no. 5, pp. 1432–1441, 2009.

J. A. Silva-Macedo, M. A. Romero, and B.-H. V. Borges, “An extended FDTD method for the analysis of electromagnetic field rotations and cloaking devices,” Progress In Electromagnetics Research, vol. 87, pp. 183–196, 2008.

Y. Luo, H. Chen, J. Zhang, L. Ran, and J. A. Kong, “Design and analytical full-wave validation of the invisibility cloaks, concentrators, and field rotators created with a general class of transformations,” Physical Review B, vol. 77, p. 125127, 2008.

W. Cai, U. K. Chettiar, A. V. Kildishev, and V. M. Shalaev, “Optical cloaking with metamaterials,” Nature Photonics, vol. 1, no. 4, pp. 224–227, 2007.

W. X. Jiang, T. J. Cui, X. M. Yang, Q. Cheng, R. Liu, and D. R. Smith, “Invisibility cloak without singularity,” Applied Physics Letters, vol. 93, no. 19, p. 194102, 2008.

G. R. Werner and J. R. Cary, “A stable FDTD algorithm for non-diagonal, anisotropic dielectrics,” Journal of Computational Physics, vol. 226, no. 1, pp. 1085–1101, 2007.

T. J. Cui, D. R. Smith, and R. Liu, Metamaterials: Theory, Design, and Applications, Springer, New York, NY, USA, 2010.

T. Allen and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Boston, MA, USA, 2005.

C. Argyropoulos, E. Kallos, and Y. Hao, “Bandwidth evaluation of dispersive transformation electromagnetics based devices,” Applied Physics A, vol. 103, no. 3, pp. 715–719, 2011.

Y. Zhao, P. Belov, and Y. Hao, “Accurate modelling of the optical properties of left-handed media using a finite-difference time-domain method,” Physical Review E, vol. 75, Article ID 037602, 2007.

Y. Yuan, K. Zhang, X. Ding, B. Ratni, S. N. Burokur, and Q. Wu, “Complementary transmissive ultra-thin meta-deflectors for broadband polarization-independent refractions in the microwave region,” Photonics Research, vol. 7, no. 1, pp. 80–88, 2019.

K. Zhang, Y. Yuan, X. Ding, B. Ratni, S. N. Burokur, and Q. Wu, “High-efficiency metalenses with switchable functionalities in microwave region,” ACS Applied Materials & Interfaces, vol. 11, no. 31, pp. 28423–28430, 2019.