Accurate and Numerically Stable FDTD Modeling of Human Skin Tissues in THz Band

We propose finite-difference time-domain (FDTD) modeling suitable for healthy skin, basal cell carcinoma, dysplastic pigmentary nevus, and non-dysplastic pigmentary nevus in the frequency range of 0.25 THz to 1.05 THz. Toward this purpose, we utilize the complex-conjugate pole-residue (CCPR) dispersion model, because it is very simple to extract the accurate CCPR coefficients using the powerful vector fitting tool. In the FDTD method, it is of great importance to check the numerical stability conditions. If the coefficients extracted through the vector fitting tool do not satisfy the numerical stability conditions, the particle swarm optimization (PSO) algorithm is employed to obtain the accurate and numerically stable coefficients. Numerical examples are provided to validate our proposed FDTD modeling.


I. INTRODUCTION
In the case of basal cell carcinoma, which accounts for 75% of all skin cancers, the incidence in Europe is increasing by 5.5% per year in recent years [1], [2]. In addition, malignant melanoma, which evolves from dysplastic pigmentary nevus, accounts for 3% of skin cancers, and has gradually increased since the 1960s, with nearly 96,000 new cases in 2019 [3]. Early detection of these skin cancers is important because they have high metastasis and poor prognosis. Terahertz (THz) imaging technology for diagnosing skin cancer has been actively studied. Note that THz radiation is non-hazardous since it is non-ionizing and low average power for THz imaging is usually used [4]. THz imaging technology of skin cancer can differentiate healthy skin tissues and diseased skin tissues [5], [6]. For example, a THz pulsed system was used to image basal cell carcinoma ex vivo for total 21 samples [7]. When basal cell carcinoma and normal tissues are compared, basal cell carcinoma has much higher absorption at THz. Therefore, the regions of basal cell carcinoma marked by THz pulsed imaging systems have a high correlation with histology. In addition, THz radiation can be harnessed for cancer treatment because of its interaction with cellular components [8], [9]. Numerical methods are The associate editor coordinating the review of this manuscript and approving it for publication was Lei Zhao . highly required to accurately analyze electromagnetic wave interaction with various human tissues in THz band.
Among many methods in computational electromagnetics [10]- [15], the finite-difference time-domain (FDTD) method has been widely used to analyze various electromagnetic wave problems due to its simplicity, robustness, and accuracy [16]- [26]. In addition, a wideband response can be obtained with a single time domain analysis. The FDTD method has been successfully used to analyze complex dispersive medium. Various dispersion models such as Debye, Drude, Lorentz, complex-conjugate pole-residue (CCPR), quadratic complex rational function (QCRF), and modified Lorentz (mLor) [27]- [38] have been introduced so far to examine electromagnetic wave interaction with dispersive media. The Debye dispersion model has been widely used for the dispersion modeling of human tissues and it was successfully used for healthy skin and basal cell carcinoma tissues in THz band [39]- [41]. In this work, the CCPR dispersion model is employed because it is a more general dispersion model of the Debye, Lorentz, and Drude models. Note that the Debye model requires one term per one pole, where both pole and residue must be real-valued and the CCPR model requires two terms (conjugate pairs) per one pole. However, the CCPR model has a higher degree of freedom than the Debye model because both pole and residue can have nonzero real and imaginary parts simultaneously, which indicating that fewer poles are required for the CCPR model than the Debye model to have the same modeling accuracy. In addition, the sum of CCPR two terms cancels out the imaginary parts due to the complex conjugate property and thus FDTD computational operation can be reduced. Therefore, the computational efficiency of CCPR-based FDTD simulations is significantly enhanced than Debye-based FDTD simulations owing to the higher degree of freedom and complex-conjugate property [33], [36]- [38]. In addition, the CCPR parameters of dispersive media can be simply obtained by the powerful vector fitting tool [42]. Based on the CCPR dispersion model, we present an accurate FDTD modeling of four different types of human skin tissues such as healthy skin [39]- [41], [43], [44], basal cell carcinoma [39]- [41], [44], dysplastic pigmentary nevus [43], [44], and non-dysplastic pigmentary nevus [43], [44] in the THz range of 0.25 THz to 1.05 THz.
Due to the explicit solution of the time-dependent Maxwell curl equations, the FDTD simulation is conditionally stable [45], [46]. Therefore, it is very important to investigate the numerical stability conditions of dispersive FDTD modeling developed in this work. It must be checked whether the extracted CCPR coefficients satisfy the numerical stability conditions of the CCPR-FDTD formulation [47]. When the extracted CCPR parameters are not numerically stable, we employ the particle swarm optimization (PSO) algorithm [48], [49] to derive the CCPR parameters that satisfy the numerical stability conditions. This additional procedure can guarantee accurate and numerically stable FDTD modeling of human skin tissues. Numerical examples are employed to validate our proposed FDTD modeling for the four human skin tissues in the frequency of interest.

II. METHODOLOGY
Assuming an e jωt time dependence, the relative permittivity of the CCPR dispersion model [33] can be expressed as where ε r,∞ is the relative permittivity at the infinite frequency, p q (rad/s) and r q (rad/s) denote the pole and residue of the CCPR model, respectively, M is the total number of the CCPR pairs, and * means the complex conjugate. As alluded previously, the aim of this work is to develop accurate and stable FDTD dispersion modeling to perform electromagnetic wave analysis on various human skin tissues. First, FDTD dispersion modeling is based on the measurement data. When the measurement data are expressed as refractive index n(ω) and absorption coefficient α(ω), they should be converted to the real part and the imaginary part of the relative permittivity using the following relations: where K is the extinction coefficient, defined as [50] K (ω) = α(ω)c 0 2ω (4) and c 0 is the speed of light in free space, c 0 = 1/ √ µ 0 ε 0 . Next, the CCPR coefficients are extracted using the simple and powerful vector fitting tool [42]. The vector fitting tool can robustly provide rational function approximations for the measured values in the frequency domain because it does not fail for poorly selected initial guess. In addition, it is very easy to implement in a computer program with the aid of standard software packages which are used to solve matrix problems.
Note that even if the dynamic stability condition is satisfied, the numerical stability conditions may not be satisfied, since the discretization in both space and time domains is not considered in the dynamic stability condition [45], [51]. The dynamic stability condition for the CCPR dispersion model is Re(p) ≤ 0, which can lead to stable and causal ε r (ω) [33]. After utilizing this vector fitting tool, in this work, we check whether the derived CCPR parameters satisfy the numerical stability conditions [47]: In above, t is the FDTD time step size, α is the FDTD space step size, and k α denotes the numerical wavenumber in the α direction [47]. In the numerical stability analysis, sin 2 ( k α α /2) is set as one for all possible numerical wavenumbers. Note that the above four numerical stability conditions are derived by the von Neumann method with the Routh-Hurwitz (R-H) criterion [27] and details can be found in [47]. The first numerical stability condition can be inferred from the dynamic stability condition and the last numerical stability condition is exactly same as the standard Courant-Friedrichs-Lewy (CFL) condition [52], [53]. In this work, if the CCPR parameters do not satisfy the above numerical stability conditions, the PSO algorithm [48], [49] is performed to extract the numerically stable coefficients. The PSO algorithm is one of the most popular optimization algorithms, which was developed with ideas from social behavior patterns of bee, bird, and fish groups. The PSO algorithm uses groups of particles called swarms and these particles navigate through the search space. Because each particle knows where it has moved before, it remembers the location with the most optimal fitness (pBest) among the positions it has passed, and shares the position with the best fitness in the swarm (gBest). By repeating this procedure, an optimal solution with the best fitness can be derived. In the PSO algorithm, the next position of each particle can be found from the following equation [49]: where v i+1 n is the velocity of particle n at iteration i + 1, x i n is the current position of particle n at iteration i, and rand() is a random number between 0 and 1. Also, c 1 is the cognitive parameter and represents confidence of particle in itself, and c 2 is the social parameter and denotes confidence of particle in its neighbors. These two parameters control the overall velocity of particle. When the value of the velocity is too large, particle may pass the position of the optimal solution, and when the velocity is too small, it may not be able to sufficiently explore the solution space. Therefore, the selection of the appropriate values of c 1 and c 2 is required. In this work, the fitness is defined as the inverse of root mean square relative error (RMSRE): where N means the number of sampling frequencies, and ε r and ε r indicate the curve-fitted data and the measurement data, respectively. In the conventional PSO algorithm, an optimal solution is derived such that it leads to the best fitness. In this work, an additional process is performed to extract the CCPR coefficients that satisfy the numerical stability conditions. In specific, when the CCPR coefficients are numerically unstable, we intentionally set the fitness value as a very small number so that the particles of the PSO algorithm can escape from the numerically unstable CCPR coefficients and find the new optimal CCPR coefficients that satisfy the numerical stability conditions. Fig. 1 shows the flow chart for CCPR-FDTD dispersion modeling developed in this work.

III. NUMERICAL EXAMPLES
As mentioned previously, we consider healthy skin, basal cell carcinoma, dysplastic pigmentary nevus, and non-dysplastic pigmentary nevus. The measurement data for refractive index n(ω) and absorption coefficient α(ω) in the frequency range of 0.25 THz to 1.05 THz can be found in [44].

A. HEALTHY SKIN AND BASAL CELL CARCINOMA
First, we extract the parameters of the CCPR dispersion model for healthy skin and basal cell carcinoma by utilizing the powerful vector fitting technique only, as the conventional approach. Let us address the total CCPR pairs number (M ). As the number of poles increases, the computational efficiency of the FDTD simulation decreases. In this work, the number of CCPR terms is determined when its RMSRE is equal to or less than 3% while gradually increasing the number of CCPR terms. The resulting CCPR parameters are listed in Table 1. It is confirmed that these extracted CCPR coefficients satisfy the numerical stability conditions of (5) as long as the FDTD time step size is limited by the CFL condition. Fig. 2 shows the complex relative permittivity of the CCPR dispersion model (solid red line) and measurement data (black dots) for healthy skin tissues. It is found out that the maximum error between the measured data and the fitted data is 4.7% and the RMSRE over the whole range of interest frequencies is 1.81%. The measured data and fitted data for basal cell carcinoma tissues are shown in Fig. 3 and its maximum error and the RMSRE are 2.5% and 1.3%, respectively. The numerically stable coefficients of the CCPR model can  be accurately derived by employing only the powerful vector fitting tool for healthy skin and basal cell carcinoma.

B. DYSPLASTIC PIGMENTARY NEVUS AND NON-DYSPLASTIC PIGMENTARY NEVUS
For dysplastic pigmentary nevus and non-dysplastic pigmentary nevus, it is found out that the CCPR parameters that satisfy the numerical stability conditions of (5) cannot be obtained by using the vector fitting tool only, different from the previous two skin tissues. To illustrate this, we first extract CCPR dispersion modeling of dysplastic pigmentary nevus by utilizing only the vector fitting tool. The resulting CCPR parameters are listed in Table 2. Fig. 4 shows that the CCPR parameters agree well with the measurement data and its RMSRE is 1.78%. Albeit with good accuracy, this CCPR dispersion modeling cannot be employed for the electromagnetic analysis of dysplastic pigmentary nevus due to numerical instability, which will be shown later.
Next, the CCPR parameters are extracted by using the vector fitting tool and the PSO algorithm simultaneously. In this work, we set the maximum number of iteration as 30, and the swarm size as 8000 for the PSO algorithm. Also, the number of independent trials is set to 30, and since the best RMSRE is derived when the cognitive parameter c 1 and the social parameter c 2 are 2, respectively, these values are used in this work. As mentioned previously, when the CCPR coefficients do not satisfy the numerical stability conditions, we deliberately set the fitness as a very small value. Table 3 shows the CCPR parameters for dysplastic pigmentary nevus and non-dysplastic pigmentary nevus by employing our proposed technique. Fig. 5 shows the measurement data (black dots) and the complex permittivity fitted through the CCPR dispersion model (red solid line) for dysplastic pigmentary  nevus. The maximum error of the fitted data against the measured data for dysplastic pigmentary nevus is 4.9% and the RMSRE over the whole frequencies is 2.05%. The plot of the measurement data of non-dysplastic pigmentary nevus and its fitted complex permittivity is shown in Fig. 6. It is observed that its maximum error and the RMSRE are 4.13% and 2.17%, respectively. Now, we perform actual FDTD simulations to investigate numerical stability of CCPR-FDTD modelings. For simplicity (without loss of generality), one-dimensional (1-D) electromagnetic analysis for dysplastic pigmentary nevus is considered as shown in Fig. 7. In the FDTD simulations, the time step size t = C n z/c 0 and the Courant number C n = 0.99 [52], [53] are employed and the FDTD grid size of z = 1.1µm is used. Note that the FDTD grid size is chosen such that the point per wavelength (PPW) is larger than 100 for the smallest wavelength in the frequency range of interest. The whole computational domain consists of 2000 FDTD cells with 1000 FDTD cells of the free space and the other 1000 FDTD cells of dysplastic pigmentary nevus. The sinewave-modulated Gaussian pulse with x-polarization is excited in the free space and denoted as where f 0 is the center frequency, t 0 is the time delay, and t w is the half width of the pulse. In this work, f 0 = 0.65 THz, t 0 = 4.43 ps, and t w = 1.48 ps to consider the frequency range of interest. Also, the 10 layers of perfectly matched layer (PML) are used for the absorbing boundary condition [16]. Fig. 8(a)   shows E x at the discontinuity interface when the numerically unstable CCPR parameters (Table 2) are used. As expected from our previous observation, the FDTD simulation suffers from numerical instability, although the CCPR parameters are in good agreement with the measurement data. Fig. 8(b) shows the E x counterpart for the numerically stable CCPR parameters (Table 3) are used. Numerical stability is observed since the CCPR parameters satisfy the numerical stability conditions of (5).
Next, we investigate the reflection coefficient for numerically stable CCPR-FDTD simulations for four different human skin tissues. FDTD simulation results for four human skin tissues and the theoretical results [44] are in good agreement with each other. The RMSRE for the reflection coefficients of the four types of skin tissues is 2.21% (healthy skin), 1.73% (basal cell carcinoma), 2.1% (dysplastic pigmentary nevus), and 1.68% (non-dysplastic pigmentary nevus). The amplitude and phase of the reflection coefficients change with the frequency. The amplitude of the    reflection coefficient is inversely proportional to the frequency in 0.25 THz -0.75 THz and the phase of the reflection coefficient is generally proportional to the frequency. The electromagnetic response depends on the skin tissue type. It is clearly observed that the amplitude slope (versus the frequency) of healthy skin tissue is the steepest and its phase slope (versus the frequency) is the flattest among four human tissues. Our proposed FDTD dispersion modeling of human skin tissues can be used to understand electromagnetic phenomena of human skin tissues and provide in-depth foundation for systematic development of THz imaging technology as well as THz therapeutics. For example, the optimal system parameters such as the operating frequency and the bandwidth can be found by various numerical simulations based on our FDTD modeling. In addition, our work can be useful for studying THz tumor diagnosis methodology such as the Cole-Cole diagram [43], [44] and its principal component analysis [43].

IV. CONCLUDING REMARKS
In the previous works [39]- [41], FDTD dispersion modeling was presented for healthy skin and basal cell carcinoma tissues but there is no literature for dysplastic pigmentary nevus and non-dysplastic pigmentary nevus. In this work, we have proposed accurate and robust FDTD modeling for healthy skin, basal cell carcinoma, dysplastic pigmentary nevus, and non-dysplastic pigmentary nevus. For accurate FDTD modeling of various human skin tissues, we employ the CCPR dispersion model since accurate CCPR parameters can be simply extracted by utilizing the vector fitting tool. First, the coefficients of the CCPR dispersion model for four different skin tissues are extracted using the vector fitting tool. Then, it is checked whether the extracted coefficients satisfy the numerical stability conditions or not. If the parameters are not numerically stable, the PSO algorithm is employed to obtain the numerically stable and accurate CCPR coefficients. It should be noted that, in this work, the numerical stability conditions are merged into the PSO algorithm so that the extracted CCPR coefficients can lead to numerical stable FDTD analysis. For healthy skin and basal cell carcinoma, the vector fitting tool is enough to extract the numerically stable CCPR coefficients and the RMSREs are less than 2%. On the other hand, the vector fitting tool and the PSO algorithm should be employed simultaneously for dysplastic pigmentary nevus and non-dysplastic pigmentary nevus to obtain accurate and robust FDTD dispersion modeling. It is observed that the RMSRE for dysplastic pigmentary nevus is 2.05% and the RMSRE for non-dysplastic pigmentary nevus is 2.17%. Actual FDTD simulations are performed to validate CCPR-FDTD modeling developed in this work. Although our accurate and robust FDTD dispersion model is applied for the CCPR dispersion model of various human skin tissues, it can be easily extended to other dispersion model of general complex dispersive media in a similar fashion. In the future, more practical 3-D examples for skin cancers will be conducted, similar to breast cancer examples [26].