International Journal of Antennas and Propagation

International Journal of Antennas and Propagation / 2020 / Article

Research Article | Open Access

Volume 2020 |Article ID 8294769 |

Ki-Woong Bae, Dong-Wook Seo, Noh-Hoon Myung, "Reduction of FDTD Stair-Casing Error regarding to Metamaterials by Using High-Order Polynomial Transformation Functions", International Journal of Antennas and Propagation, vol. 2020, Article ID 8294769, 16 pages, 2020.

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

Academic Editor: Giuseppe Castaldi
Received09 Nov 2019
Accepted24 Dec 2019
Published13 Jan 2020


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 [35]. 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 [79], were designed by a linear transformation. In [79], 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 [712], 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:wherexi 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 [79] 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, R1 and R2 (R1 and R3 for the concentrator).

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 f1(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 R2 and expands the origin into R1. In mathematical expressions,  = 1 at r′ = R2 and , respectively. Thus, the conditions should be  = 1 and  = 0, where means . Next, since the rotator does not rotate fields in r′ > R2 and r′ < R1, 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 R3 and with the core at R1, which is expressed as  = 1 and  = R1/R2, 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 × R1 = 1.5 × R2 = R3. These are summarized in Table 1.



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

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 f3 are continuous, although those with the use of f1 are discontinuous.

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 R1 < r′ < R1 +  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′ = R1 +  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 f5 = f3. 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 f5, 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 f5 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 f5 and its material parameter change near boundaries are smaller than those obtained with f3. Notice that the parameter change in the middle is larger with the use of f5 than that with the use of f3, as shown in Figure 5.

Cloak ()Rotator ()Concentrator ()

A7.6312e + 04−4.4973e + 055.7501e + 04
B−5.9825e + 043.4810e + 05−7.5034e + 04
C1.8069e + 04−1.0356e + 053.8622e + 04
D−2.6118e + 031.4726e + 04−9.7941e + 03
E1.8148e + 02−1.0059e + 031.2247e + 03

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 R2 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.

R1 (mm)R2 (mm)R3 (mm)f0 (GHz)x (=∆y)


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):

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.

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.

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.

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.

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.

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 R1 is also included in the graph.

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.


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


  1. 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. View at: Publisher Site | Google Scholar
  2. 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. View at: Publisher Site | Google Scholar
  3. 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. View at: Publisher Site | Google Scholar
  4. 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. View at: Publisher Site | Google Scholar
  5. 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. View at: Publisher Site | Google Scholar
  6. 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. View at: Publisher Site | Google Scholar
  7. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science, vol. 312, no. 5781, pp. 1780–1782, 2006. View at: Publisher Site | Google Scholar
  8. H. Chen and C. T. Chan, “Transformation media that rotate electromagnetic fields,” Applied Physics Letters, vol. 90, no. 24, p. 241105, 2007. View at: Publisher Site | Google Scholar
  9. 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. View at: Publisher Site | Google Scholar
  10. 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. View at: Publisher Site | Google Scholar
  11. 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. View at: Publisher Site | Google Scholar
  12. 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. View at: Google Scholar
  13. 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. View at: Publisher Site | Google Scholar
  14. 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. View at: Publisher Site | Google Scholar
  15. 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. View at: Publisher Site | Google Scholar
  16. 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. View at: Publisher Site | Google Scholar
  17. T. J. Cui, D. R. Smith, and R. Liu, Metamaterials: Theory, Design, and Applications, Springer, New York, NY, USA, 2010.
  18. T. Allen and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Boston, MA, USA, 2005.
  19. 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. View at: Publisher Site | Google Scholar
  20. 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. View at: Publisher Site | Google Scholar
  21. 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. View at: Publisher Site | Google Scholar
  22. 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. View at: Publisher Site | Google Scholar

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.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

Article of the Year Award: Outstanding research contributions of 2020, as selected by our Chief Editors. Read the winning articles.