Perfectly Matched Layer Formulation of the INBC-FDTD Algorithm for Electromagnetic Analysis of Thin Film Materials

The impedance network boundary condition (INBC)-based finite-difference time-domain (FDTD) method has been widely used for electromagnetic analysis of highly conductive thin film materials. In the INBC-FDTD formulation, the electromagnetic field variations inside the thin film material are taken into account mathematically and thus extremely small FDTD grids are not necessary for the FDTD modeling of the material. Therefore, computational efficiency of the INBC-FDTD formulation is significantly better than other FDTD formulations. Albeit with this great advantage, the INBC-FDTD formulation cannot be fully employed for thin film materials because the corresponding perfectly matched layer (PML) formulation has not been reported in literature. In this work, we propose a PML formulation suitable for the INBC-FDTD algorithm. Numerical examples illustrate that the proposed PML-INBC-FDTD formulation can yield good absorption performance and also it can improve computational efficiency while maintaining numerical accuracy.


I. INTRODUCTION
The finite-difference time-domain (FDTD) method has been popularly employed for a variety of research areas including dispersion-engineered metamaterials, plasma, photonics, and biomedical applications [1]- [7]. The main advantages of the FDTD method are simplicity, robustness, and broadband analysis. The FDTD method can be used to analyze electromagnetic wave propagation in complex media such as nonlinear, dispersive, and chiral media [8]- [12].
Thin film materials have been widely used for electromagnetic applications [13]- [16]. When modeling a thin film material in the standard FDTD method, the electromagnetic field variations within the material must be taken into account accurately. Therefore, in the standard FDTD method, very refined spatial grids are usually used, which leads to overwhelming computational resources [17]. To tackle this problem, the impedance network boundary condition The associate editor coordinating the review of this manuscript and approving it for publication was Yi Ren .
(INBC)-FDTD formulation was presented [18]. The INBC algorithm considers the tangential electric and magnetic field components located on the surface of the thin film material and employs the impedance relationships between them. The change in electromagnetic fields inside the thin film material is considered mathematically in the INBC algorithm. Therefore, in the INBC-FDTD formulation, very refined spatial grids are not needed any more but the conventional FDTD grids can be employed, which leads to improve computational efficiency. The INBC-FDTD formulation was successfully used to analyze shielding effectiveness of highly conductive thin film materials [19]- [21]. In [19], [20], the INBC-FDTD modeling of thin conductive shields was proposed. In [21], thin composite multilayered panels were efficiently analyzed using the INBC-FDTD formulation.
Since the computational simulation space is finite, an absorbing boundary condition must be used for the outermost boundaries of the FDTD discrete space. Among various absorbing boundary conditions, the perfectly matched layer (PML) is robust and powerful because it VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ can effectively absorb arbitrarily-incident and polarized electromagnetic waves in complex media without spurious reflections [22]- [24]. However, a PML formulation for the INBC-FDTD formulation has not been reported in literature.
In this work, we propose a PML formulation suitable for the INBC-FDTD algorithm. In the next section, the INBC algorithm is briefly reviewed and then the formulation for the thin film materials is presented. Next, numerical examples are employed to illustrate that the proposed PML-INBC-FDTD formulation can effectively absorb electromagnetic waves. Finally, concluding remarks are followed.

II. PML-INBC-FDTD FORMULATION
As shown in Fig. 1, let us consider a thin film material with thickness d along the x-axis in y = J y on the two-dimensional xy-plane for simplicity without loss of generality. Here, ξ indicates the FDTD grid size along the ξ direction. The tangential electromagnetic field components on the surface of the thin film material are E x and H z . According to the INBC algorithm [19], the tangential magnetic field on the upper (or lower) surface is expressed as H + z (or H − z ). Similarly, the tangential E + x and E − x components are introduced on the surface. In the INBC algorithm, two impedance relationships are considered to connect tangential electromagnetic field components on both sides of the thin film material by the self and mutual impedances: where Z 11 and Z 22 indicate the self impedance, and Z 12 and Z 21 represent the mutual impedance. The self impedance and the mutual impedance are given by where η = µ ε + σ/jω , Note that ε, µ, and σ indicate the permittivity, the permeability, and the conductivity of the thin film material respectively. In above, the intrinsic impedance is denoted as η, and the propagation constant is expressed as γ . Note that the thickness (d) of the thin film material is taken into account mathematically in (2). As alluded previously, the INBC-based FDTD method can use the conventional spatial grid size, which can lead to dramatically increase its computational efficiency [19]. In the FDTD method, it is necessary to convert the frequency-domain impedance relationships into the time-domain counterparts. However, it is not straightforward to derive the INBC-FDTD formulation by directly using the inverse Fourier transform (IFT) of (2). Alternatively, in the INBC-FDTD algorithm, (2) is approximated as the series of the rational functions by utilizing the vector fitting method proposed in [25]: Note that (3) can be simply converted into the time-domain functions via the IFT.
As mentioned previously, the PML absorbing boundary condition for the INBC-FDTD formulation has not been presented until now. The PML absorbing boundary condition was applied to only the conventional FDTD formulation except the INBC-FDTD formulation for the thin film material [19]. Note that this conventional INBC-FDTD formulation can lead to acceptable numerical results for thin copper film in 10 MHz-10 GHz [19] but it cannot be employed for some EM analyses of thin film materials due to spurious PML reflections.
In this work, we propose the PML formulation suitable for the INBC-FDTD algorithm in order to fully employ the INBC-FDTD algorithm for electromagnetic analysis of thin film materials. The complex-frequency-shifted (CFS)-PML is employed in this work and the CFS-PML stretching variable is written as [26] where κ ξ increases attenuation of evanescent waves, and α ξ ensures late-time stability for low frequency waves. The PML can be implemented through modified Maxwell's equations with stretched coordinates [27], [28]. In the INBC-FDTD formulation, only Ampere's law is considered for tangential electromagnetic fields [19] and thus the modifed Ampere's law (in the frequency domain) in the PML region can be written as Inserting (4) into the above equations, we can write where g − zx , g + zx , g − zy , and g + zy are auxiliary PML variables and they are defined by By applying the IFT to the above equations, the auxiliary differential equations can be written as The FDTD update equations for H ± z can be derived by using the IFT and the finite difference schemes to (6): Here, t denotes the FDTD time step size, the subscript indicates the space index, and the superscript indicates the time index. The grid position of each field component is shown in Fig. 1. Note that the PML auxiliary variables g − zx and g − zy are collocated with H − z and the PML auxiliary variables g + zx and g + zy are collocated with H + z . According to the INBC-FDTD algorithm [19], the forward or backward difference scheme to the normal directional spatial derivative (e.g., the y-directional spatial derivative in Fig. 1). The backward difference scheme (BDS) is applied to ∂ ∂y E − x in (9) and the forward difference scheme (FDS) is utilized to ∂ ∂y E + x in (9). Following the CFS-PML implementation described in [9], [26], the update equations for the PML auxiliary variables can be written as where τ = κ ξ ε 0 κ ξ α ξ + σ ξ .
It should be noted that we employ the BDS in (13) and the FDS in (14) for the consistency with the INBC-FDTD algorithm.

III. NUMERICAL EXAMPLES
Numerical examples are employed to validate the proposed PML-INBC-FDTD algorithm. As a first example, we consider graphene, a carbon allotrope that forms a two-dimensional planar structure in the form of a hexagonal grid. Graphene is attracting attention as a key material in various research areas such as photodetectors, optical sensings, and metamaterials [29]- [31]. Since graphene can yield tunable surface plasmon polariton (SPP) in the THz band, it is suitable for compact THz electromagnetic wave devices [32]- [35]. The surface conductivity (σ s ) of graphene is described by the Kubo's formula [36]. When the thickness of graphene is d, the volume conductivity is obtained by σ v (ω) = σ s (ω)/d and thus the relative complex permittivity can be written as ε r (ω) = 1 + σ s (ω)/jωε 0 d [10], [37]. In this work, the thickness of graphene is set to 0.33 nm [38].
By using the vector fitting method, we extract four pairs of pole-residue values in (3) for graphene in the frequency range of 1-30 THz [25]. Table 1 lists the extracted values of poles and residues for the self impedance of graphene and Table 2 lists the extracted values of poles and residues for the mutual impedance of graphene. The self impedance and mutual impedance obtained by the vector fitting method agree very well with the analytical impedances, as shown in Fig. 2 and Fig. 3.
The relative reflection errors of the conventional INBC-FDTD and the proposed PML-INBC-FDTD formulations are calculated to compare absorption performance. We consider two parallel graphene sheets with the separation of 50 nm, as shown in Fig. 4. The FDTD grid size is 1 nm, the time step is set to 1 × 10 −18 s for the Courant-Friedrichs-Levy (CFL) stability condition, and a Gaussian-modulated sinewave with the bandwidth of 1-30 THz band is used for source excitation [10], [37]. The whole FDTD domain consists of 120 × 120 cells including 10-cell PML layers in each direction. The y-directed electric current sources are excited and then E y at the corner of the two parallel graphene sheets is observed. To obtain a reference solution, a very large FDTD domain with 64-cell PML layers is used [39]. All FDTD simulations are performed on Intel i7-10700 CPU. As shown in Fig. 5, the proposed PML-INBC-FDTD formulation can yield good absorption performance while the absorption performance of the conventional INBC-FDTD formulation is poor. Note that this graphene parallel plate waveguide can yield SPP wave propagation and thus electromagnetic fields on the graphene sheets within the PML region are not negligible.
It is also interesting to investigate field distribution obtained by both INBC-FDTD simulations. For this study, TABLE 1. Extracted values of poles(p 11,k = p 22,k = p s,k ) and residues(r 11,k = r 22,k = r s,k ) for the self impedance of graphene where z 11,∞ = 297.5853.

TABLE 2.
Extracted values of poles (p 12,k = p 21,k = p m,k ) and residues(r 12,k = r 21,k = r m,k ) for the mutual impedance of graphene where z 12,∞ = 297.5737.    setup is maintained as the previous example. The electric field at the time step of 2.4 × 10 5 for the conventional INBC-FDTD simulation is illustrated in Fig. 6. As shown in the figure, the electric field is severely contaminated by spurious reflections. Fig. 7 shows the electric field at the same time step for the proposed PML-INBC-FDTD simulation. It is observed that there is no spurious reflection on the computational boundaries, implying again that the proposed PML-INBC-FDTD formulation works very well.
Let us consider the accuracy of the proposed PML-INBC-FDTD formulation quantitatively. Toward this purpose, we calculate the propagation constant of the graphene parallel plate waveguide by using the discrete Fourier transform (DFT) of the FDTD time data along the waveguide, addressed in [10]. Fig. 8 shows the propagation constant extracted by the proposed PML-INBC-FDTD simulation and the theoretical propagation constant [38]. The proposed FDTD result is consistent with the theoretical counterpart. Next, we employ larger FDTD grid sizes, i.e., 2 nm, 5 nm, and 10 nm while the thickness of graphene is maintained  as 0.33 nm. As shown in Fig. 8, all PML-INBC-FDTD simulations are in good agreement with the theoretical values. As alluded previously, the change in electromagnetic fields inside graphene is naturally considered in the PML-INBC-FDTD formulation and thus very refined spatial grids are not necessary. It should worthy noting that very refined FDTD grid size (e.g., 1 nm) should be employed for other PML-FDTD formulations without utilizing the INBC algorithm [10], [37]. Table 3 summarizes CPU time and memory storage for the PML-INBC-FDTD simulation with four different FDTD grid sizes. As the FDTD grid size increases, computational efficiency increases dramatically while maintaining computational accuracy. The proposed PML-INBC-FDTD simulation with the FDTD grid size of 10 nm needs 46 times less CPU time and costs 23 times less memory storage, compared to the 1-nm FDTD case.
As the next example, we consider 400 µm-thick carbon fiber composites (CFC) with the relative permittivity of 6.4 and conductivity of 1.5 × 10 4 [40] in the frequency range of 100 Hz-10 GHz and the geometry of this problem is shown in Fig. 9. The FDTD grid size is 1 mm and the CFL number of 0.99 is used. The conventional INBC-FDTD simulation and the proposed INBC-FDTD simulation are   performed in 70 × 70 FDTD cells including 10-cell PML layers in each direction. The z−directed magnetic current with a differentiated Gaussian pulse in the time domain is excited and then H z is observed. In addition, the FDTD simulation with the extremely refined FDTD grid size, i.e., s = 10 µm, is performed for the same geometry to obtain the reference solution. The normalized magnetic fields of the three numerical results are illustrated in Fig. 10. The proposed PML-INBC-FDTD result is in good agreement with