An Efficient 3-D Stochastic HIE-FDTD Algorithm for Investigation of Statistical Variation in Electromagnetic Field

We propose a stochastic hybrid implicit–explicit finite-difference time-domain method (S-HIE-FDTD) to compute the mean and variance of the electromagnetic (EM) fields using a single simulation, given those of the conductivity and permittivity in the computation domain. The mean and variance field update equations underlying the proposed method are derived from the field update equations of the “traditional” deterministic HIE-FDTD. The Courant–Friedrichs–Lewy condition of the S-HIE-FDTD depends on the spatial discretization sizes only in two dimensions; therefore, for computation domains with fine geometric features only in the remaining dimension, it uses a time step size that is larger than that of fully explicit schemes. Indeed, numerical results demonstrate that the proposed method is faster than the previously developed stochastic FDTD in computing the mean and variance of the EM fields in two different problems: wave propagation through a multilayer human tissue and transmission through a frequency selective surface.


I. INTRODUCTION
The statistical electromagnetic (EM) analysis [1] is indispensable in various fields of engineering and physics, including, but not limited to, bioelectromagnetics [2]- [4], atmospheric wave propagation [5], remote sensing [6], [7], and integrated circuit design [8], where materials' electrical properties are not deterministically known (often due to measurement limitations) and/or scatterer and antenna dimensions are only known within a range of values (often due to fabrication tolerances). The purpose of such an analysis is to predict the mean and variance [and, sometimes, the probability distribution function (pdf)] of a quantity of interest (QoI) (e.g., induced voltage and transmitted power) when the mean and variance (and the pdf) of input/simulation parameters (e.g., conductivity, permittivity, scatterer, and antenna dimensions) are provided.
A common way to statistically characterize a QoI is the Monte Carlo (MC) method [9]. The MC method samples the domain of input parameters and computes the QoI at every sampling point using a deterministic EM solver. These samples of the QoI are then used to compute its mean, variance, and pdf. The MC method is easy to implement since it does not call for any modifications to existing EM solvers (i.e., it is nonintrusive), but it requires a large number of samples (EM solver executions) to produce accurate results.
To avoid sampling in the domain of input parameters, recently, a novel stochastic finite-difference time-domain method (S-FDTD) has been developed [10]. This method is built upon the deterministic FDTD [11]- [13], but it is reformulated and implemented to produce the mean and variance of the EM fields given those of the conductivity and permittivity in the simulation domain. However, just like its deterministic counterpart, the time step size of the S-FDTD is determined by the finest spatial discretization size through the Courant-Friedrich-Levy (CFL) condition [11], [14] but not the desired level of accuracy, especially for structures that have fine and coarse geometric features at the same time.
To this end, various unconditionally stable (deterministic) FDTD schemes [15]- [17] have been developed to remove the CFL restriction. These schemes can use larger time step sizes, but their computational cost per time step is increased due to the implicit field updates. Therefore, they might unnecessarily increase the overall computational cost for problems involving structures with fine geometric features only in one or two dimensions (e.g., radiation from patch antennas, transmission through thin layers, and scattering in a layered medium). To address this shortcoming, several FDTD schemes, such as weakly conditionally stable (WCS) method [18], [19] and hybrid implicit-explicit method [20]- [23], have been developed. These schemes use implicit updates only in one or two dimensions (ideally where fine geometric features exist). This has two consequences: 1) the spatial discretization size along these dimensions does not "contribute" to the CFL condition, which results in a larger time step size compared with fully explicit schemes and 2) their computational cost per time step is lower compared with fully implicit schemes. Therefore, the overall computational cost of these schemes is expected to be lower than that of the fully explicit or fully implicit FDTD methods when analyzing structures with fine geometric features only in one or two dimensions.
In this communication, we develop a stochastic hybrid implicitexplicit FDTD (S-HIE-FDTD) to efficiently compute the mean and variance of the EM fields given those of the conductivity and permittivity in a simulation with fine geometric features in one or This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ two dimensions. To come up with the update equations of the S-HIE-FDTD, we reformulate the update equations of the deterministic HIE-FDTD to take the mean and variance of the conductivity and permittivity as inputs. A short stability analysis shows that the S-HIE-FDTD has a CFL condition similar to that of the HIE-FDTD, and as a result, the S-HIE-FDTD has the same computational benefits as its deterministic counterpart as described in the previous paragraph. Indeed, the numerical results demonstrate that the S-HIE-FDTD is more efficient than the S-FDTD in the statistical analysis of wave propagation through a multilayer biological tissue and transmission through a frequency selective surface (FSS).

II. FORMULATION
The full set of equations, which are solved/updated by the HIE-FDTD for all electric and magnetic field components, can be found in [23]. In this communication, only the equations, which update the x-component of the electric field and the z-component of the magnetic field, namely, E x and H z , are used to demonstrate how the HIE-FDTD is adopted to compute the mean and variance of the fields given those of the relative permittivity ε r and conductivity σ . The HIE-FDTD updates E x and H z using Here, ε 0 and μ 0 are the permittivity and permeability in free space, ν 1 = (2ε 0 ε r − σ t)/(2ε 0 ε r + σ t), ν 2 = 2 t/(2ε 0 ε r + σ t), ν 3 = t/μ 0 , t is the time step, and x, y, and z are the space steps. In (1) and (2) and in the rest of the text, the superscript n and the subscript i, j, k mean that the pertinent field variables are sampled at time n t and location (i x, j y, k z).

A. Stochastic Model
Let g(X 1 , . . . , X N ) represent a function of random variables X 1 , . . . , X N with mean values μ X 1 , . . . , μ X N . Using the linearity of where a i are constants, the Delta method [24], and following the derivation in [10], we obtain: The variance of g(X 1 , . . . , X N ) is expressed as where S[.] and S 2 [.] represent the standard deviation and variance operators. Similarly, using the Delta method [24] and following the derivation in [10], we obtain: In addition, the following identity is available for the variance operation [10]: is the covariance and ρ X i ,X j is the correlation coefficient between X i and X j .
Update (1) and (2), and equalities (3)-(6) are used together to derive the update equations for the S-HIE-FDTD, namely, the mean and variance field equations, as described in Sections II-B and II-C, respectively. This derivation assumes that the samples of E x , E y , H y , and H z , namely, i+1/2, j,k+1/2 , and H z | n i+1/2, j +1/2,k , the permittivity ε r , and the conductivity σ are random variables. In addition, it is assumed that the same type (electric or magnetic) field samples separated by a single time step and/or a single space step in any direction (including those multiplied with coefficients ν 1 , ν 2 , or ν 3 ) are highly correlated, i.e., the relevant correlation coefficient is equal to 1 [10] ρ As a result, the mean and variance field update equations of the proposed S-HIE-FDTD involve only the correlation coefficients between field samples, ε r , and σ (i.e., ρ ε r ,E , ρ σ,E , ρ ε r ,H , ρ σ,H , and ρ ε r ,σ ). The values of these coefficients can be predicted using the methods described in [25]- [27].

C. Variance Field Equations
Let the function g(·) be defined by (1). The first step in the derivation of the variance field equation for E x is to take the variance of both sides of (1) Applying (6) to the left-hand side of (10) and using (7) while noticing ν 1 is a function of random variables ε r and σ , we obtain Applying (5) and using (7), we obtain (5) and (6) multiple times to the right-hand side of (10) and using (7), we obtain In (12) and (13) where s ε r and s σ are the standard deviations of ε r and σ . Adding μ 2 ν 2 κ 2 2 E[h] 2 to the right-hand side of (13) and subtracting it yield Assuming (14) and taking the square-root of the resulting equation, we obtain Using (11) and (12) in the left-hand side of (10) and (15) in its right-hand side, we obtain In (16), the standard deviation values of the field samples are the unknowns to be updated during marching in time.
We follow a similar procedure to obtain the variance field equation for H z . Let the function g(·) be defined by (2). The first step is to take the variance of both sides of (2) Applying (6) to the left-and right-hand sides of (17) and using (7) in the resulting equation, we obtain Note that ν 3 is a constant (not a function of a random variable).
Taking the square root of (18) and rearranging the terms yield In (19), the standard deviation values of the field samples are the unknowns to be updated during marching in time.

D. Numerical Stability Analysis
For the sake of simplicity of the numerical stability analysis, the conductivity is assumed deterministic and set to zero. Following the Fourier method [10], the CFL condition of the S-HIE-FDTD is where c denotes the speed of light in free space. This CFL condition uses only the mean value of the permittivity but not its variance.
In other words, the CFL of the S-HIE-FDTD is the same as that of the deterministic HIE-FDTD with permittivity replaced with its mean.

III. NUMERICAL EXAMPLES
In this section, the accuracy and efficiency of the S-HIE-FDTD are compared with those of a traditional MC method [9] and/or the S-FDTD [10] for two numerical examples. In both the examples, absorbing boundaries enforce the first-order standard Mur absorbing boundary condition (ABC), as described in [10] and [11], and the plane wave excitation is introduced using the total-field/scatteredfield (TF/SF) technique described in [11]. All the computations are performed on a desktop computer with a 32-GB RAM and a 3.90-GHz Intel Core i7-3770 processor.

A. Wave-Propagation Through Multilayer Biological Tissue
For the first example, we use the proposed method to analyze EM wave propagation through human tissue consisting of three layers: skin, fat, and muscle (see Fig. 1). The mean and standard deviation of the conductivity and permittivity of these layers are given in Table I. The simulation domain used in this problem is shown in Fig. 1. The dimensions of the simulation domain are x ∈ [0, 10] mm, y ∈ [0, 500] mm, and z ∈ [0, 10] mm. The ABCs and the periodic boundary conditions (PBCs) are used along the y-, x-, and z-directions, respectively. A plane wave propagating along the y-direction with x-polarized and unit-amplitude electric field is incident on the tissue layers. The time signature of the excitation is a sine function with a frequency of 2 GHz. Two sets of simulations are carried out using this setup. For both of the sets, the S-HIE-FDTD and S-FDTD are executed for ρ = 0.5 and ρ = 1. Here, ρ is the correlation coefficient between the permittivity/conductivity and the field samples.
In the first set, all three layers of tissue have the same thickness of 54 mm. The spatial mesh uses x = y = z = 1 mm and t = 6.67 × 10 −4 ns for the S-HIE-FDTD and S-FDTD. The MC method executes the (deterministic) FDTD for 10 000 times. For each execution, the permittivity and conductivity of the tissue layers are selected independently of the Gaussian distributions with the mean and standard deviation, as provided in Table I. The S-HIE-FDTD, S-FDTD, and MC method are used to compute the mean and variance of the electric field's x-component at x = 5 mm, y ∈ [0, 500] mm, and z = 5 mm. Fig. 2(a) and (b) shows the  plots of the amplitude and phase of this mean (after transformed into frequency domain at 2 GHz) versus y and demonstrates that all three methods produce practically the same mean values everywhere in the simulation domain. The fact that the result from the MC method matches those obtained by the S-FDTD and S-HIE-FDTD shows that the assumption of the high correlation between the samples of the electric field and between the samples of the magnetic field [see (7)] is accurate. Fig. 3 shows the plots of the time-domain amplitude of the variance versus y. The results obtained by the S-HIE-FDTD and S-FDTD agree very well with each other for both values of ρ. However, results with ρ = 0.5 agree better with the result obtained by the MC method. We believe this means that setting ρ = 1 overestimates the correlation.
In the second set of simulations, the thickness of the skin, fat, and muscle layers is 5.4, 54, and 42 mm, respectively. The spatial mesh uses x/3 = y = z/3 = 0.5 mm and t = 3 × 10 −3 ns for the S-HIE-FDTD and t = 1 × 10 −3 ns for the S-FDTD. Note that y is smaller than x and z because the layers are located along the y-direction and skin layer is very thin. Also, note that the CFL condition of the S-HIE-FDTD has no constraint in the y-direction, that is why t of the S-HIE-FDTD does not have to be reduced and is larger than t of the S-FDTD. The rest of the simulation parameters is the same as those used in the first set of simulations.  Table II presents the computational requirements of the MC method, S-HIE-FDTD, and S-FDTD. The S-FDTD and S-HIE-FDTD are significantly faster than the MC method. Here, the number of (total) iterations is equal to the number of time steps multiplied by the number of executions. The number of executions is 1 for the S-FDTD and S-HIE-FDTD and 10 000 for the MC method. The S-HIE-FDTD is roughly 1.7 times faster than the S-FDTD with only a 20% increase in the memory requirement. The memory requirement of the MC method is significantly larger than those of the S-HIE-FDTD and S-FDTD since it stores field samples needed to compute the mean and variance.
We note here that the Fourier transform of the variance of a timedomain field is not equal to the variance of its frequency-domain counterpart (variance is nonlinear in field values). That is why we plot the variance in the time domain in Figs. 3 and 5. This is also the approach adopted in several other communications [10], [25].

B. Transmission Through an FSS
In the second example, we use the proposed method to analyze transmission through an FSS [21] supported by a dielectric substrate (see Fig. 6). The top view of FSS unit cell is shown in Fig. 6(a),   where D = 14 mm, g = 10 mm, a = 10 mm, and w = 1 mm, and the thickness of the substrate is 2 mm. The simulation domain used in this problem is shown in Fig. 6(b). The dimensions of the simulation domain are x ∈ [0, 160] mm, y ∈ [0, 10] mm, and z ∈ [0, 14] mm. The ABCs and PBCs are used along the x-, y-, and z-directions, respectively. A plane wave propagating along the x-direction with ypolarized and unit-amplitude electric field is incident on the FSS. The time signature of the excitation is given by G(t) = e −4π(t −t 0 ) 2 /τ 2 , where τ = 6.67 × 10 −2 ns and t 0 = 0.8τ .
The spatial mesh uses x/5 = y = z/5 = 0.2 mm and t = 2.22 × 10 −3 ns for the S-HIE-FDTD and t = 5.56×10 −4 ns for the S-FDTD. Note that y has to be smaller than x and z since the slot on the surface of the FSS is very thin along the y-direction. Also, note that the CFL condition of the S-HIE-FDTD has no constraint in the y-direction, that is why t of the S-HIE-FDTD is larger than that of the S-FDTD.
Four sets of simulations are carried out, where S-FDTD and S-HIE-FDTD are executed for four different values of the substrate permittivity's mean μ ε r ∈ {1.9, 2.2, 2.4, 2.6}. Its variance is set to 0.1 for all simulations. The mean and the standard deviation of the transmission coefficients are computed in the frequency domain from the results obtained using the S-FDTD and S-HIE-FDTD. Let "M" and "S" denote the mean and the standard deviation, respectively. Fig. 7 shows the plots of M, M-S, and M+S obtained using the S-FDTD and S-HIE-FDTD (a total of six curves) with μ ε r = 2.2. The results obtained by the two methods match very well. Moreover, the figure shows that the uncertainty in the permittivity  of the substrates results in a shift in the resonance frequency of the transmittance. Fig. 8 shows the plots of M obtained using the S-HIE-FDTD with μ ε r ∈ {1.9, 2.2, 2.4, 2.6}. A red/blue shift is observed with the increase/decrease of the permittivity that coincides with the results of M+S/M−S. It should be noted that the shift of the peak in Fig. 8 is more obvious than that in Fig. 7. The reason is that μ ε r is fixed, and the plus and minus S (see Fig. 7) only leads to a weak influence compared with the change of mean (see Fig. 8). Table III presents the computational requirements of the S-HIE-FDTD and S-FDTD and shows that the S-HIE-FDTD is roughly 2.3 times faster with only 20% increase in the memory requirement.

IV. CONCLUSION
An S-HIE-FDTD to compute the mean and variance of the EM fields by a single simulation given those of the conductivity and permittivity in the computation domain is presented. The mean and variance field update equations are derived. The CFL condition of the S-HIE-FDTD method resembles that of the HIE-FDTD method, resulting in a time step size larger than that would be used by an explicit scheme.
Two numerical examples are provided to show the accuracy, efficiency, and applicability of the proposed S-HIE-FDTD. In the first example, EM wave propagation through a human tissue consisting of three layers is analyzed. In each layer, the permittivity and conductivity are defined by their mean and standard deviation. In the second example, transmission through an FSS supported by a dielectric substrate is analyzed. The permittivity of this substrate is defined by its mean and variance. The results show that the proposed method can dramatically reduce the computation time without any reduction in the accuracy level compared with the S-FDTD method. However, the variance of the EM fields still demonstrates a difference with respect to results obtained using the MC method; we believe this is due to the approximation of the correlation coefficient between the field samples and the permittivity/conductivity.