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

Authors: Ki-Woong Bae; Dong-Wook Seo; Noh-Hoon Myung
Journal: International Journal of Antennas and Propagation
 (2020)

---

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


---

## Body

## 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 = μ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 Section2. 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:(1)εi′j′=detgi′j′−1/2gi′j′ε,(2)μi′j′=detgi′j′−1/2gi′j′μ,where(3)gi′j′=Λki′Λlj′δkl,(4)Λii′=∂xi′∂xi=∂fxi∂xi,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.Iff ′ is continuous at a boundary with the derivative of space before the transformation ∂xi/∂xi, that is, Λii′ = I, then parameter matching is achieved as εi′j′ = ε and μi′j′ = μ. 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):(5)g1r′=r=r′−R1R2−R1R2,R1≤r′≤R2,(6)f1r=−1R2−R1r+R2R2−R1,R1≤r≤R2,where g means the inverse function of f. These linear transformation functions are continuous at both boundaries, R1 and R2 (R1 and R3 for the concentrator).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.
(a)
(b)
(c)(7)
f
1
r
=
R
1
R
2
r
,
0
≤
r
≤
R
2
,
R
3
−
R
1
R
3
−
R
2
r
−
R
2
−
R
1
R
3
−
R
2
R
3
,
R
2
<
r
≤
R
3
,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 functionf1(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:(8)f3r=r′=Ar3+Br2+Cr+D,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, ∂r′/∂r = 1 at r′ = R2 and limr′⟶R1∂r′/∂r⟶∞, respectively. Thus, the conditions should be g3′R2 = 1 and g3′R1 = 0, where g3′ means ∂g3r′/∂r′. Next, since the rotator does not rotate fields in r′ > R2 and r′ < R1, the rotation rate f3′ should be 0 at both boundaries, which can be written as f3′R1 = f3′R2 = 0. Finally, the concentrator is in contact with the background at R3 and with the core at R1, which is expressed as f3′R3 = 1 and f3′R1 = R1/R2, respectively. These conditions are, respectively, summarized as follows:(9)g3R1=0,g3R2=R2,g3′R1=0,g3′R2=1,(10)f3R1=1,f3R2=0,f3′R1=0,f3′R2=0,(11)f3R2=R1,f3R3=R3,f3′R2=R1R2,f3′R3=1.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:(12)R33R32R31R23R22R213R322R3103R222R210ABCD=f3R3f3R2f3′R3f3′R2=R3R11R1/R2.The unknown coefficients were derived with the radii given as 3 ×R1 = 1.5 × R2 = R3. 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.73918The graphs of linear and third-order functions are compared in Figure2. For example, the linear function f1 in Figure 2(a) has discontinuous derivatives, such as f1′R1=f1′R2=0.5, while f3 shows continuous derivatives, such as g3′R1=1/∞⟶0 and g3′R2=1. Note that f3 is monotonic in range 0 < r < R2; thus, g3′=1/f3′.Figure 2
Transformation functions: (a) the cloak, (b) the rotator, and (c) the concentrator.
(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 f3 are continuous, although those with the use of f1 are discontinuous.Figure 3
Metamaterial parameters: (a) the cloak, (b) the rotator, and (c) the concentrator.
(a)
(b)
(c)


### 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 rangesR1 < r′ < R1 + 2Δx and R2,3 −2Δx < r′ < R2,3 has the possibility of being situated in partially filled cells. Thus, it is logical to set conditions at r′ = R1 + 2Δx and R2,3 − 2Δx, where Δx denotes the FDTD spatial step size. Six (that is, 4 + 2) conditions can be met by using a fifth-order polynomial function:(13)f5r=Ar5+Br4+Cr3+Dr2+Er+F.Iff5′R1+2Δx = f3′R1+2Δx and f5′R2,3−2Δx = f3′R2,3−2Δx, then f5 = f3. If f5′R1+2Δx = f3′R1 and f5′R2,3−2Δx = f3′R2,3, then the change is completely suppressed because f5′R1 = f3′R1 and f5′R2,3 = f3′R2,3. Thus, it is logical to choose f5′R1+2Δx value in the range of f3′R1+2Δx and f3′R1. In (11), 0 ≤ P ≤ 1 is used to express f5′R1+2Δx and f5′R2,3−2Δx conveniently. Note that the change is more suppressed with increasing P:(14a)f5′R1+2Δx=f3′R1+2Δx+f3′R1−f3′R1+2Δx×P,(14b)f5′R2,3−2Δx=f3′R2,3−2Δx+f3′R2,3−f3′R2,3−2Δx×P.The functionf5, on the other hand, is nonmonotonic with high P values. This can be explained by Rolle’s theorem and inflection points. A detailed explanation will be omitted for brevity. The upper limits of P that make f5 monotonic were calculated by increasing P from 0.01 to 0.99. The calculation reveals that the maximum P values are 0.70, 0.55, and 0.63, for the cloak, the rotator, and the concentrator, respectively.CoefficientsA through F were derived from a 6 × 6 matrix similar to (12) with the aforementioned P 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.Table 2
Coefficients of fifth-order polynomial functions.
Cloak (P=0.7)
Rotator (P=0.55)
Concentrator (P=0.63)
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.4071Figure 4
Proposed fifth-order transformation function: (a) the cloak, (b) the rotator, and (c) the concentrator. The function had smaller changes near boundaries compared withf3.
(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 withf3.
(a)
(b)
(c)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. Figure6 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 εr,εΦ,andμz of the three MTMs in consideration. The smoother parameter change indicates that it will produce smaller FDTD stair-casing errors.Figure 6
Theεr component of the field rotator around the R2 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).
(a)
(b)
(c)Table 3
FDTD simulation conditions.
R
1 (mm)
R
2 (mm)
R
3 (mm)
f
0 (GHz)
∆x (=∆y)
Cloak
106
212
—
2
λ
0/30
Rotator
106
212
—
2
λ
0/30
Concentrator
106
212
318
2
λ
0/30The proposed high-order function method is applicable to elliptic coordinate systems, as shown in Figure7(a) [15]. For example, consider a cloak in an elliptic system with a focal length p=0.08m. 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):(15)ξ′=fξ=Aξ3+Bξ2+Cξ+D.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.
(a)
(b)
(c)The unknown coefficients were derived from the continuity conditions off and f ′ as follows: A = –8.47518, B = 44.70059, C = –74.56999, and D = 40.05591.The material parameters are provided in Figure7(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.


## 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:(8)f3r=r′=Ar3+Br2+Cr+D,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, ∂r′/∂r = 1 at r′ = R2 and limr′⟶R1∂r′/∂r⟶∞, respectively. Thus, the conditions should be g3′R2 = 1 and g3′R1 = 0, where g3′ means ∂g3r′/∂r′. Next, since the rotator does not rotate fields in r′ > R2 and r′ < R1, the rotation rate f3′ should be 0 at both boundaries, which can be written as f3′R1 = f3′R2 = 0. Finally, the concentrator is in contact with the background at R3 and with the core at R1, which is expressed as f3′R3 = 1 and f3′R1 = R1/R2, respectively. These conditions are, respectively, summarized as follows:(9)g3R1=0,g3R2=R2,g3′R1=0,g3′R2=1,(10)f3R1=1,f3R2=0,f3′R1=0,f3′R2=0,(11)f3R2=R1,f3R3=R3,f3′R2=R1R2,f3′R3=1.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:(12)R33R32R31R23R22R213R322R3103R222R210ABCD=f3R3f3R2f3′R3f3′R2=R3R11R1/R2.The unknown coefficients were derived with the radii given as 3 ×R1 = 1.5 × R2 = R3. 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.73918The graphs of linear and third-order functions are compared in Figure2. For example, the linear function f1 in Figure 2(a) has discontinuous derivatives, such as f1′R1=f1′R2=0.5, while f3 shows continuous derivatives, such as g3′R1=1/∞⟶0 and g3′R2=1. Note that f3 is monotonic in range 0 < r < R2; thus, g3′=1/f3′.Figure 2
Transformation functions: (a) the cloak, (b) the rotator, and (c) the concentrator.
(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 f3 are continuous, although those with the use of f1 are discontinuous.Figure 3
Metamaterial parameters: (a) the cloak, (b) the rotator, and (c) the concentrator.
(a)
(b)
(c)


## 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 rangesR1 < r′ < R1 + 2Δx and R2,3 −2Δx < r′ < R2,3 has the possibility of being situated in partially filled cells. Thus, it is logical to set conditions at r′ = R1 + 2Δx and R2,3 − 2Δx, where Δx denotes the FDTD spatial step size. Six (that is, 4 + 2) conditions can be met by using a fifth-order polynomial function:(13)f5r=Ar5+Br4+Cr3+Dr2+Er+F.Iff5′R1+2Δx = f3′R1+2Δx and f5′R2,3−2Δx = f3′R2,3−2Δx, then f5 = f3. If f5′R1+2Δx = f3′R1 and f5′R2,3−2Δx = f3′R2,3, then the change is completely suppressed because f5′R1 = f3′R1 and f5′R2,3 = f3′R2,3. Thus, it is logical to choose f5′R1+2Δx value in the range of f3′R1+2Δx and f3′R1. In (11), 0 ≤ P ≤ 1 is used to express f5′R1+2Δx and f5′R2,3−2Δx conveniently. Note that the change is more suppressed with increasing P:(14a)f5′R1+2Δx=f3′R1+2Δx+f3′R1−f3′R1+2Δx×P,(14b)f5′R2,3−2Δx=f3′R2,3−2Δx+f3′R2,3−f3′R2,3−2Δx×P.The functionf5, on the other hand, is nonmonotonic with high P values. This can be explained by Rolle’s theorem and inflection points. A detailed explanation will be omitted for brevity. The upper limits of P that make f5 monotonic were calculated by increasing P from 0.01 to 0.99. The calculation reveals that the maximum P values are 0.70, 0.55, and 0.63, for the cloak, the rotator, and the concentrator, respectively.CoefficientsA through F were derived from a 6 × 6 matrix similar to (12) with the aforementioned P 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.Table 2
Coefficients of fifth-order polynomial functions.
Cloak (P=0.7)
Rotator (P=0.55)
Concentrator (P=0.63)
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.4071Figure 4
Proposed fifth-order transformation function: (a) the cloak, (b) the rotator, and (c) the concentrator. The function had smaller changes near boundaries compared withf3.
(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 withf3.
(a)
(b)
(c)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. Figure6 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 εr,εΦ,andμz of the three MTMs in consideration. The smoother parameter change indicates that it will produce smaller FDTD stair-casing errors.Figure 6
Theεr component of the field rotator around the R2 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).
(a)
(b)
(c)Table 3
FDTD simulation conditions.
R
1 (mm)
R
2 (mm)
R
3 (mm)
f
0 (GHz)
∆x (=∆y)
Cloak
106
212
—
2
λ
0/30
Rotator
106
212
—
2
λ
0/30
Concentrator
106
212
318
2
λ
0/30The proposed high-order function method is applicable to elliptic coordinate systems, as shown in Figure7(a) [15]. For example, consider a cloak in an elliptic system with a focal length p=0.08m. 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):(15)ξ′=fξ=Aξ3+Bξ2+Cξ+D.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.
(a)
(b)
(c)The unknown coefficients were derived from the continuity conditions off and f ′ as follows: A = –8.47518, B = 44.70059, C = –74.56999, and D = 40.05591.The material parameters are provided in Figure7(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 Figure8 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 Figures9 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.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.
(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:(16)Relative RMSE=1/K∑FDTD cellKHzanalytic−HzFDTD2∑FDTD cellKHzanalytic2,where K is the number of FDTD cells. The error was proportional to O1/Δx2 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.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.
(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 Figure11(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 Figure12. 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.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.
(a)
(b)
(c)As seen in Figure5, 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 Figure1(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.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.
(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. TheHz 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.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.
(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 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.



---

*Source: 8294769-2020-01-13.xml*