Research Article  Open Access
Reduction of FDTD StairCasing Error regarding to Metamaterials by Using HighOrder Polynomial Transformation Functions
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 highorder polynomial functions as the transformation function. Since similar materials are filled in boundary cells of the finitedifference timedomain (FDTD) algorithm, the staircasing 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 finitedifference timedomain (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 staircasing approximation. The partially filled cells require special treatments [1, 2], such as subpixel smoothing [3–5]. The staircasing approximation fills a boundary cell with one of its filling media; the staircasing 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, staircase 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 secondorder 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. Highorder 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 staircasing 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:wherex^{i} and x^{i′} 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 ∂x^{i}/∂x^{i}, that is, = I, then parameter matching is achieved as = ε and = μ. Here, f ′ is the derivative of f with respect to x^{i} 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_{1} and R_{2} (R_{1} and R_{3} for the concentrator).
(a)
(b)
(c)
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_{1}(r) = Ar + B has two coefficients, A and B; thus, the four conditions cannot always be met. On the other hand, an Nthorder polynomial function has N + 1 coefficients. These coefficients can be freely adjusted to meet the required conditions.
3. Proposal of HighOrder Polynomial Functions
3.1. ThirdOrder Polynomial Functions
The aforementioned four conditions can be met by using a thirdorder 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_{2} and expands the origin into R_{1}. In mathematical expressions, = 1 at r′ = R_{2} and , respectively. Thus, the conditions should be = 1 and = 0, where means . Next, since the rotator does not rotate fields in r′ > R_{2} and r′ < R_{1}, 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_{3} and with the core at R_{1}, which is expressed as = 1 and = R_{1}/R_{2}, 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} = 1.5 × R_{2} = R_{3}. These are summarized in Table 1.

The graphs of linear and thirdorder functions are compared in Figure 2. For example, the linear function f_{1} in Figure 2(a) has discontinuous derivatives, such as , while f_{3} shows continuous derivatives, such as and . Note that f_{3} is monotonic in range 0 < r < R_{2}; thus, .
(a)
(b)
(c)
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_{3} are continuous, although those with the use of f_{1} are discontinuous.
(a)
(b)
(c)
3.2. FifthOrder Polynomial Functions
The FDTD staircasing 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 staircasing 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_{1} < r′ < R_{1} + and R_{2,3} − < r′ < R_{2,3} has the possibility of being situated in partially filled cells. Thus, it is logical to set conditions at r′ = R_{1} + and R_{2,3} − , where Δx denotes the FDTD spatial step size. Six (that is, 4 + 2) conditions can be met by using a fifthorder polynomial function:
If = and = , then f_{5} = f_{3}. 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_{5}, 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_{5} 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_{5} and its material parameter change near boundaries are smaller than those obtained with f_{3}. Notice that the parameter change in the middle is larger with the use of f_{5} than that with the use of f_{3,} as shown in Figure 5.

(a)
(b)
(c)
(a)
(b)
(c)
The staircasing error in the FDTD algorithm originates from the abrupt material parameter change. Hence, exactly how the staircasing 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 staircased R_{2} 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 thirdorder function in Figure 6(b). The fifthorder 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 staircasing errors.
(a)
(b)
(c)

The proposed highorder 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 thirdorder function similar to (8):
(a)
(b)
(c)
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 staircased 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 thirdorder 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 secondorder FDTD method and truncated with the secondorder Mur absorbing boundary condition. The MTMs were staircasing 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 thirdorder 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 staircasing 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 fifthorder function, which is attributable to the smaller material parameter change near boundaries.
(a)
(b)
(c)
(a)
(b)
(c)
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 secondorder FDTD algorithm was used [18]. The fifthorder function had the smallest error followed by the thirdorder and linear functions; the highorder 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)
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 fifthorder function had larger time increase than the thirdorder 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)
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 highorder 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 thirdorder 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)
(a)
(b)
(c)
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 H_{z} field was spatially averaged along a line parallel to the yaxis. 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_{1} is also included in the graph.
(a)
(b)
(c)
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 higherorder 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 thirdorder and fifthorder functions. The higher transmitted amplitude in Figure 15(a) indicates that the higherorder 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 staircasing 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 highorder 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 highorder polynomials. As a result, the FDTD staircasing 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 higherorder polynomials that have smaller material change near boundaries. These results imply that the proposed method would permit accurate FDTD analysis of complexshaped 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, “Contourpath effective permittivities for the twodimensional finitedifference timedomain method,” Optics Express, vol. 13, no. 25, pp. 10367–10381, 2005. View at: Publisher Site  Google Scholar
 J. Nadobny, D. Sullivan, and P. Wust, “A general threedimensional tensor FDTDformulation for electrically inhomogeneous lossy media using the Ztransform,” IEEE Transactions on Antennas and Propagation, vol. 56, pp. 1027–1040, 2008. View at: Publisher Site  Google Scholar
 A. F. Oskooi, C. Kottke, and S. G. Johnson, “Accurate finitedifference timedomain simulation of anisotropic media by subpixel smoothing,” Optics Letters, vol. 34, no. 18, pp. 2778–2780, 2009. View at: Publisher Site  Google Scholar
 A. Deinega and I. Valuev, “Subpixel smoothing for conductive and dispersive media in the finitedifference timedomain method,” Optics Letters, vol. 32, no. 23, pp. 3429–3431, 2007. View at: Publisher Site  Google Scholar
 J. Liu, M. Brio, and J. V. Moloney, “Subpixel smoothing finitedifference timedomain 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
 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
 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
 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
 W. X. Jiang et al., “Design of arbitrarily shaped concentrators based on conformally optical transformation of nonuniform rational Bspline surfaces,” Applied Physics Letters, vol. 92, no. 26, p. 264101, 2008. View at: Publisher Site  Google Scholar
 Y. Cui, C. Argyropoulos, and Y. Hao, “Fullwave finitedifference timedomain simulation of electromagnetic cloaking structures,” Optics Express, vol. 16, no. 9, pp. 6717–6730, 2008. View at: Publisher Site  Google Scholar
 C. Argyropoulos, Y. Zhao, and Y. Hao, “A radiallydependent dispersive finitedifference timedomain 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
 J. A. SilvaMacedo, 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
 Y. Luo, H. Chen, J. Zhang, L. Ran, and J. A. Kong, “Design and analytical fullwave 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
 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
 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
 G. R. Werner and J. R. Cary, “A stable FDTD algorithm for nondiagonal, anisotropic dielectrics,” Journal of Computational Physics, vol. 226, no. 1, pp. 1085–1101, 2007. View at: Publisher Site  Google Scholar
 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 FiniteDifference TimeDomain 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. View at: Publisher Site  Google Scholar
 Y. Zhao, P. Belov, and Y. Hao, “Accurate modelling of the optical properties of lefthanded media using a finitedifference timedomain method,” Physical Review E, vol. 75, Article ID 037602, 2007. View at: Publisher Site  Google Scholar
 Y. Yuan, K. Zhang, X. Ding, B. Ratni, S. N. Burokur, and Q. Wu, “Complementary transmissive ultrathin metadeflectors for broadband polarizationindependent refractions in the microwave region,” Photonics Research, vol. 7, no. 1, pp. 80–88, 2019. View at: Publisher Site  Google Scholar
 K. Zhang, Y. Yuan, X. Ding, B. Ratni, S. N. Burokur, and Q. Wu, “Highefficiency 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
Copyright © 2020 KiWoong 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.