A Short Review of FDTD-Based Methods for Uncertainty Quantification in Computational Electromagnetics

We provide a review of selected computational methodologies that are based on the deterministic finite-difference time-domain algorithm and are suitable for the investigation of electromagnetic problems involving uncertainties. As it will become apparent, several alternatives capable of performing uncertainty quantification in a variety of cases exist, each one exhibiting different qualities and ranges of applicability, which we intend to point out here. Given the numerous available approaches, the purpose of this paper is to clarify the main strengths and weaknesses of the described methodologies and help the potential readers to safely select the most suitable approach for their problem under consideration.


Introduction
Uncertainty quantification commonly emerges as a significant and multiaspect issue in numerous engineering and other scientific disciplines.In the context of electromagnetic (EM) problems, randomness may appear in several forms, affecting, for example, geometric features, electric parameters, and input sources.Even though it is a rather common approach to treat EM problems as completely deterministic (and several reliable time-and frequency-domain computational algorithms have been developed for their solution), there exist several instances where one or more parameters/features/aspects of an EM problem cannot be strictly determined.This lack of knowledge apparently induces a degree of randomness in the output quantities, which in many cases can be extremely important.Scientists and engineers have realized the significance of being able to reliably describe a system's or problem's variability and have proposed a number of different treatments for such cases.The scope of this review paper is to examine some of the available solutions, in the context of computational electromagnetics, providing information for both standard and more modern approaches.Evidently, the pertinent research area is quite vast, and a complete review of all available techniques in the context of a single work is quite difficult, if not impossible.Considering that, in the area of time-dependent EM problems, perhaps the most popular simulation method is the finite-difference timedomain (FDTD) algorithm [1], we decided to turn our attention to stochastic methods that rely on, or are combined with, the FDTD approach.Most of the cited works in this paper are quite contemporary, and we believe that the following analysis will clarify the strengths and weaknesses of the considered approaches, as well as the level of suitability for specific types of applications.

Monte Carlo Methods
Monte Carlo (MC) methodologies [2] are commonly used for the study of stochastic EM problems, as they provide a standard line of work that, in most cases, can lead to reliable results.The main advantage of MC approaches is the fact that they can be implemented easily, without usually requiring exceptional skills or advanced mathematical knowledge.In fact, the algorithmic development of MC-based solutions is rather straightforward and, certainly, less complicated than other existing approaches.On the other hand, solutions via MC techniques normally entail augmented computational cost (e.g., due to repeated implementation of full-wave numerical methods) and converge rather slowly.Nevertheless, thanks to their reliability, MC solutions are often used as a means to assess the validity of other approaches and can be easily applied to diverse scientific disciplines, without requiring significant modifications.In the following paragraphs, we present a small fraction of examples utilizing MC methods for the investigation of time-dependent EM problems via FDTD schemes.
Back to 1991, the FDTD method was applied in the context of MC simulations for the study of EM scattering from random rough surfaces [3].Specifically, considering a Gaussian surface with specific root-mean-square height and correlation length, the average scattered intensity at different scattering angles was obtained, after performing 50 FDTD runs.This low number of samples should be probably attributed to the limited computing resources available at the time.The results were compared against data from frequency-domain approaches, and good agreement was observed.The wellknown advantage of using a time-domain method to simultaneously obtain results at different frequencies was also noted.In a similar context, the work in [4] examined the wave scattering from surfaces with either Gaussian or multiscale Pierson-Moskowitz roughness spectra, using standard as well as contour-path FDTD methods, to deal with the unwanted staircase phenomenon.The performance of these two methods was evaluated by computing the radar cross section for different scattering angles, and the detrimental effect of staircasing was found to be more prominent at higher values of incident angles.
In [5], the 3D FDTD method was combined with MC analysis for the investigation of scattering from perfectly conducting objects, when the latter are embedded in random inhomogeneous soils with spatially varying permittivity.The MC analysis was specifically used to obtain the bistatic radar cross section for normal plane-wave incidence, considering a total of 100 runs.Coherent and incoherent averagings were performed to extract the scattering of the object alone and the contributions of the medium fluctuations, respectively.The scattering from inhomogeneous media was also discussed in [6], where the two-dimensional FDTD method and the MC approach were combined.Scattering patterns were deduced, after extracting the coherent and incoherent (diffused) parts of the scattered energy.A homogenization procedure was also proposed, and comparison between lossy and lossless cases was performed.
The combined MC-FDTD method was applied in [7] for the extraction of the frequency-dependent complex conductivity of doped silicon, and much more accurate results compared to existing Drude models were derived.A completely different problem was solved in [8], where the electricfield intensity was estimated close to human body models in realistic environments.Various exposure scenarios and different frequencies were considered, and a total of 50,000 simulations were conducted for each case.In [9], the distribution of light induced by an electron beam in a fluorescent thin film was estimated via the MC-FDTD scheme.In this case, uncertainty pertains to the determination of the step length that an electron travels between successive scatterings in order to compute the electron trajectory.Finally, a complete study of uncertainties characterizing lightning-generated EM fields was performed in [10].The three-dimensional curvilinear FDTD method [1] was used to investigate the effects of ground inhomogeneities, terrain roughness, geometric irregularity of the lightning channel, and channel's inclination.Due to the high computational burden involved, 60-100 simulations were conducted in each case, and computations were accelerated by up to 100x via code parallelization on graphics processing units [11].

Stochastic Finite-Difference Time-Domain Algorithms
The methods of this section are variants of the standard FDTD algorithm [1], with similar structure and computational complexity.It is reminded that the implementation of the FDTD method requires a Cartesian grid, where electricand magnetic-field components comply with a staggered spatial arrangement.A similar staggering approach is applied in time as well.This spatiotemporal arrangement enables the efficient implementation of second-order finite-difference approximations for the first-order Maxwell's system, which produce an explicitly solved set of difference equations.For instance, the update of the   component in a lossless medium with ,  constitutive parameters is performed via where |  ,, stands for (Δ, Δ, Δ, Δ), Δ denotes the selected discretization step with respect to the  variable ( ∈ {, , , }), ,  correspond to electric-and magnetic-field components, and  describes source current-density components.Five more equations similar to the aforementioned one complete the system of FDTD update equations in threedimensional EM problems.

The Stochastic FDTD Method.
A variant of the standard FDTD algorithm which is capable of predicting statistical moments, such as the mean value E{⋅} and the standard deviation {⋅} of the EM fields due to uncertainty in the electric properties of the material involved in a specific setup, is developed in [12].The so-called Stochastic FDTD algorithm preserves much of the structure and features (e.g., Cartesian grid and explicitness) of the conventional deterministic approach and directly provides update equations for the aforementioned moments of the involved components.In essence, considering a function  that depends on a number of random variables  1 ,  2 , . . .,   , a second-order approximation of the mean value of  is obtained simply by E { ( 1 ,  2 , . . .,   )} ≃  ( 1 ,  2 , . . .,   ) , where  1 ,  2 , . . .,   are the mean values of  1 ,  2 , . . .,   .This result implies that the update equations of the deterministic algorithm can be used for the calculation of the expected values of the field components, without other modifications, provided that the mean electric properties of the background materials (that are the sources of uncertainty) are used in the update equations.
The extraction of the update equations for the standard deviation is trickier, as it involves applying  2 {⋅} on both sides of the deterministic finite-difference equations.Now, standard formulae such as must be used, where   is the correlation coefficient of , .At a first step, this coefficient pertains to successive (in time) values of the same field components; hence it can be safely set to 1.However, in the case of products of random quantities, a different line of work must be followed.Specifically, the following formula must be eventually applied: Unavoidably, now the correlation coefficient between field values and variables describing the material uncertainty at the same position appears.Its exact value is unknown, and we usually select (as a simple solution) again  = 1, although this choice has the tendency to overestimate the computational outcome, as it practically corresponds to a (probably unrealistic) worst-case scenario.The final form of the update equations for the standard deviation is derived after applying the square root to both sides and taking a secondorder Taylor approximation.For reference, the update of the standard deviation of   in the case of one-dimensional problems becomes where   ,   ,   , and   are constants depending, among others, on the statistical properties of the involved medium.Evidently, one needs to update the expected values of the EM components first and then proceed to the calculation of the standard deviations at each point of the grid.This technique has been applied in [12] to a scattering problem concerning a three-layer biological tissue.Although the outcome was not exact, when compared to reference data from MC, the authors showed that they could provide reliable bounds for the solution by performing two Stochastic FDTD simulations with different choices of the correlation coefficients.

Stochastic FDTD Method for Magnetized Plasma.
As the original Stochastic FDTD method was first presented in onedimensional formulation, a three-dimensional extension is developed in [13].Specifically, this version of the stochastic algorithm is suitable for problems involving wave propagation in anisotropic magnetized plasma, which can be found in cases of signal propagation within the Earth's ionosphere.Unlike [12], the method of [13] deals with Maxwell's equations coupled to current equations derived from the Lorentz equation of motion.The source of uncertainty in this case is the electron density, which exhibits a highly complex behavior in ionosphere.The same principles and approximations as the ones already described for the original Stochastic FDTD algorithm are applied.Again, it is directly found that the original difference equations can be used for the average of the variables.In fact, an overall of 18 equations is deduced.The first subset of 6 equations involves the magnetic-field components and their deviations, while the second subset comprises the components of the electric-field intensity and the current density.The last group of equations involves the deviations of E and J.The tests that were conducted considered the case of collisionless plasma, and the comparisons against MC references verified that the mean value of fields can be calculated with good accuracy, while the reliability of the standard deviation strongly depends on the selected value of the (unknown) correlation coefficient.As in the case of the original stochastic method, selecting  = 1 overestimated the results, while a much smaller value (0.05 for the considered problem) improved the accuracy significantly.

Curvilinear Stochastic FDTD Method.
Another threedimensional extension of the Stochastic FDTD algorithm is developed in [14], where a generalized curvilinear formulation is presented.The starting point now is the curvilinear FDTD method [1], which involves the update equations for the covariant and contravariant components, as well as the necessary metrics-weighted interpolation processes.The final update equations are derived by applying the methodology already described.For instance, an update equation for the deviation of the first covariant electric-field component becomes where covariant components are identified by subscripts and contravariant components by superscripts and   ,   ,   ,   ,   , and   are constants depending on material parameters and may involve metric coefficients.The method is also adapted to model geometric uncertainties by properly taking into account the statistics of metric coefficients.Finally, computations are accelerated with the aid of parallel programming on graphics processing units [11].The curvilinear Stochastic FDTD technique is applied therein to problems with uncertainty in the material electric properties, geometric length, or geometric curvature.The accuracy of the results is in agreement with the conclusions of the already mentioned works.The same approach has been already used for the solution of diverse problems including the investigation of lightning-generated EM fields over nondeterministic terrains [15] and nanomaterials and graphene-based configurations [16].Especially in the last case, a combination of different spatial stencils leads to an optimum finite-difference approximation for spatial derivatives, while the incorporation of a domain-decomposition approach results in improved versatility.

Unconditionally Stable Stochastic FDTD Methods. As
(almost) all explicit FDTD schemes, the aforementioned methodologies are conditionally stable; that is, the allowable time-step size is bounded.This restriction can result in excessive computational times when simulating problems involving wave propagation over long distances or highly dense grids, which are necessary for modeling structures with fine geometric features.Unconditionally stable methods based on special temporal updating, such as the locally-onedimensional (LOD) approach, are not restricted by stability requirements and allow the utilization of any time-step size (of course, the larger the time step, the higher the numerical errors).A LOD-based Stochastic FDTD method is presented in [17], which is suitable for the solution of the telegraph equations for transmission-line problems.Similar to the deterministic algorithm, the updates for the mean values and standard deviations of the involved variables (voltages and currents) are implicit and require the solution of tridiagonal systems at each time step.However, the additional computational burden is compensated for by the reduced number of iterations, thanks to the larger values of time steps.The numerical results in [17] verify that the statistics of wave phenomena along transmission lines that display uncertainty in their per-unit parameters can be successfully assessed even with large time step, provided that a sufficient spatial resolution is ensured.

A Possible Improvement.
It has become apparent that a nontrivial limitation of the Stochastic FDTD algorithm emanates from the lack of knowledge of certain correlation coefficients.Until today, no efficient means of their calculation has been proposed, rendering the selection of random values (commonly equal to 1) the most usual choice.A solution that partially sacrifices the simplicity and rather low complexity of the original Stochastic FDTD approach is developed in [18], where a hybrid scheme that incorporates a fraction of MC simulations attempts to remedy the aforementioned situation.Considering that correlation coefficients are calculated from where Cov stands for the covariance operator, [18] proposes their approximate estimation after performing a small number of runs, which should be enough to guarantee a relatively reliable output.The authors therein realize that the correlation coefficients also exhibit temporal not just spatial variation, which apparently poses additional difficulties.However, for time-harmonic problems, they claim that it suffices to use constant-valued coefficients, especially those calculated at the last iteration of the algorithm, when fields have finally converged.For the 1D simulations performed in [18], it appears that the suggested line of work is successful, even when the correlation coefficients are extracted from only 100 MC simulations (which correspond to only 1% of the runs required to obtain reference data).These concepts are also extended to two-dimensional problems in [19], where the variance of the specific absorption rate in the human head due to plane-wave illumination is computed.Again, considerable reduction of the error of the Stochastic FDTD method is verified (although it should be mentioned that the actual problem setup is not two-dimensional in reality).After testing different choices of the required MC runs, it was found that satisfactory convergence of the correlation coefficient is accomplished with only 25 simulations, as no significant changes were noted when performing additional tests.

An Alternative Stochastic, Single Realization, FDTD
Scheme.A stochastic realization of the FDTD algorithm, which is suitable for problems involving weak scattering phenomena, is described in [20].The authors therein realized that the (unknown) correlation coefficients appearing in the updates of the Stochastic FDTD method are due to the utilization of Taylor approximations for multiplicative-noise terms (i.e., random variables related through multiplication).
Instead, an iterative procedure is proposed, which turns the multiplicative noise into additive noise that is easier to handle.Furthermore, the proposed approach can be implemented as a single-realization algorithm, unlike MC solutions.This method was shown to perform well in weak scattering problems, where the impacts of the selected correlation lengths and the spectral content of the incident waves on the accuracy of the predicted mean field values were investigated as well.

Methods Based on Polynomial Chaos
Polynomial Chaos (PC) enables the representation of random processes with the aid of convergent polynomial series [21].
If, for example,  depends on the  independent random variables  1 ,  2 , . . .,   , then its PC expansion takes the form where  ℓ are constant coefficients,  = [ 1  2 ⋅ ⋅ ⋅   ] T is the vector of the random variables, and Ψ ℓ are proper basis functions.In practice, the aforementioned series is truncated to include  + 1 terms, and an approximation of the exact quantity is obtained.If  is the polynomial order of the expansion, then the number of terms is determined by +1 = (+)!/(!!).The basis functions are constructed as follows: that is, as products of univariate polynomials ( ℓ  is the multiindex that corresponds to the expansion order [22], ℓ = 0, . . ., ).The basis functions should be selected in accordance with the type of variable that they describe: Gaussian variables require Hermite polynomials, Jacobi polynomials are suitable for beta distributions, uniform variables require Legendre polynomials, and so forth.Note that the basis functions are orthogonal: where the inner product is defined as () is the joint probability-density function, and Γ  is the dimensional random space.After the PC approximation of a random quantity has been obtained, the corresponding mean value and variance are directly computed from the available expansion coefficients: There exist two major categories of uncertainty quantification algorithms that are based on PC expansions: intrusive and nonintrusive.In the first category, an existing (deterministic) algorithm is modified, so that it directly computes the necessary expansion coefficients.Nonintrusive schemes exploit deterministic methods to achieve their goal and require multiple executions of these methods, according to specific samplings of the random space.Although nonintrusive solutions are simpler to implement, as no changes to existing algorithms are required, intrusive approaches have the advantage of providing the expansion coefficients in the context of a single run.

Intrusive PC-FDTD Methods.
The PC-FDTD method [23] is a FDTD algorithm whose update equations directly calculate the expansion coefficients of the PC series rather than the actual fields.The pertinent equations can be derived easily from the difference equations of the original scheme after substituting each field component by the corresponding PC expansion and then applying Galerkin projection on every basis function.Taking into account the orthogonality property of the basis functions, a different update equation is obtained for each expansion coefficient.For instance, if we assume that implying that the uppercase letters denote the random field components, while lowercase letters are used for the expansion coefficients, then ( 1) is transformed to which determines the -th expansion coefficient, when  exhibits uncertainty.If, on the other hand, the corresponding node does not involve a material with variability, then the update equations of the expansion coefficients are less complex and similar to the original ones.Evidently, the computational cost of the PC-FDTD method is higher than that of the deterministic algorithm, as it involves  + 1 updates for each field component.Due to the well-known issue of "curse of dimensionality," the number of coefficients increases significantly when the number of random variables is considerable.On the other hand, especially for low or moderate number of random variables, the PC-FDTD method is far more efficient than MC approaches, which normally require thousands of simulations to provide convergent results.
In the original work of [23], the PC-FDTD method is applied to an EM compatibility problem and to plane-wave scattering from a dielectric sphere.Some later contributions described more advanced applications of the considered approach.For example, the purpose of [24] is to assess the uncertainty in wireless indoor channel models due to the randomness in the electric parameters of the wall materials of an office building or the geometric dimensions of the walls themselves.The necessary integrals for the calculation of the inner products were computed using sparse Smolyak grids [25] and Gauss-Legendre integration rules to improve efficiency.Several interesting results were extracted, including probability-density functions of the path loss along selected sectors, relative contribution of each material property to the overall loss at different spatial locations, and spatial density plots of 95% confidence intervals.
Another application of the PC-FDTD method is described in [26], where the impact of geometric uncertainties of microwave circuits on their performance is investigated.Now, the spatial steps of the finite-difference discretization are treated as random variables; however, the resulting PC-FDTD equations are very similar to, for example, (14).Many cases are treated by simply distorting the cell size of the rectilinear Yee grid; however, more general uncertainties, pertinent to statistically independent variables, call for local mesh distortions.To deal with the latter case, a curvilinear PC-FDTD method is developed, which requires that the inner products are different for each (distorted) cell.On the other hand, outside the distorted curvilinear part of the mesh, the  + 1 update equations are decoupled, as in the standard PC-FDTD technique.The proposed method was applied to various configurations, including a low-pass microstrip filter, a cascaded stub line filter, and a directional coupler, and statistics as well as sensitivities were extracted for timeand frequency-domain quantities.
The PC-FDTD method has been also applied to problems involving complex wave propagation in the ionosphere [27], similar to the Stochastic FDTD algorithm of [13].Given that only one random variable was considered in that problem (the electron density), up to fourth-order PC expansions were tested, without increasing the computational cost prohibitively.It was deduced that second-order expansions are sufficient for mean-value calculations, whereas standard deviations require at least fourth-order expansions to capture all the necessary details.Finally, the effect of ground material uncertainties on lightning-induced EM pulses was the subject of [28], where the necessary polynomial orders for sufficiently reliable outcomes were investigated, and results regarding the computation of probability-density functions and Sobol indices were presented.Due to the symmetry of the considered problem, a rotationally symmetric twodimensional FDTD algorithm was used as the starting point for the development of the stochastic approach.

Nonintrusive PC Methods.
As already mentioned, the main disadvantage of the methods of the previous subsection pertains to the modifications that a deterministic algorithm has to undergo in order to be adapted to stochastic problems.Alternatively, one can exploit the deterministic solver in a different manner [29] for the computation of the expansion coefficients involved in (8).For instance, spectral projection enables the direct calculation of the unknowns from which stems from the orthogonality of the basis functions.Evidently, the -dimensional integral needs to be computed numerically, following some rule, for example, where  () are selected integration points and  () are the corresponding weights.There exist various possibilities regarding the aforementioned numerical integrations, which include full tensor-product grids (easily obtained from onedimensional rules), sparse grids (e.g., Smolyak [25]), and MC formulae [30].It is noted that full grids are characterized by an exponentially growing total number of nodes, while the sparse solutions involve a not as fast increasing computational cost.
where  1 ,  2 , . . .,   are the corresponding deterministic outputs.It is of crucial importance that the number of equations is not smaller than the number of unknowns; hence the number of grid nodes should satisfy  ≥  + 1.The overdetermined system can be solved approximately via the leastsquares method; that is, Evidently, the aforementioned line of work is the same regardless of the selected deterministic solver.Hence, in the rest of this section, we mention a few representative EM applications, where the FDTD method is selected as the deterministic solver.For example, the variability of the specific absorption rate in human body models is examined in [31], when factors such as the plane-wave incidence angle or the position of a mobile wireless device exhibit uncertainty.A completely different subject is discussed in [32], which focuses on the characterization of wireless channels in burning buildings.In this case, flames are modeled as cold plasma; hence a dispersive FDTD algorithm is required.The Kronrod-Patterson quadrature is used for the calculation of the necessary inner products, and the uncertainty characterizing the recorded power at different positions is assessed, when there is randomness in the plasma parameters.In [33], the influence of the ground's material parameters on the EM pulses produced by lightning flashes is explored with stochastic collocation methods, utilizing either full or sparse grids.The modeling of complex stochastic problems appearing in the field or electromagnetic compatibility is discussed in [34], and problems involving the shielding effectiveness of cabinets and a dielectric resonator antenna are solved.

Discussion
Although the present review is short and cannot be considered by any means exhaustive, we believe that we have included some very representative works that represent distinct numerical approaches pertinent to uncertainty quantification in computational electromagnetics.In this concluding section, we revise the most basic features of each presented category of methods.
The major strength of MC-based approaches is their straightforward implementation, which involves direct use of existing deterministic solvers.On the other hand, their convergence is rather slow (but independent of the number of random variables) and requires a considerable amount of samples in order to provide reliable results.This becomes a nontrivial problem when full-wave solvers are used for the calculation of the deterministic solutions.Furthermore, it becomes clear that the Stochastic FDTD methods discussed in Section 3 are quite easy to implement, and the involved computational cost is just twice the cost of the corresponding deterministic algorithm.The lack of knowledge of the necessary correlation coefficients is a nontrivial deficiency, which can be partly remedied by performing a small fraction of MC simulations.Of course, such a solution increases the complexity of implementation.In addition, the Stochastic FDTD schemes calculate only the mean value and the standard deviation of random quantities.In other words, significant measures such as the probability-density functions or sensitivity indices cannot be computed, which is probably the price we have to pay for such a low-complexity approach.PC-based methodologies are more powerful in the sense that they provide a more consistent description of random processes.Intrusive PC-FDTD schemes have the advantage of calculating the PC representations in the context of a single simulation.Unfortunately, they require significant modifications of the deterministic algorithms, and the involved complexity is certainly higher.Nonintrusive methods, on the other hand, are simpler to implement, as they directly apply the deterministic scheme.However, they do need special samplings of the random space, in order to calculate the expansion coefficients reliably, in the context of either numerical quadrature or a least-squares approach.Finally, PC-based approaches suffer from the "curse of dimensionality": unlike MC methods, their efficiency deteriorates as the number of random variables increases.This means that if the dimensionality of the random space is large, PC methods do not provide any practical advantage over MC approaches in terms of the required number of simulations.In any case, each methodology exhibits certain features and advantages that should be taken into consideration, when their suitability for a certain application needs to be determined.Uncertainty quantification is a vast research field with applications in several diverse scientific disciplines, and various other approaches that are not necessarily related to the FDTD method exist, which have not been included in the present review.Some indicative examples include Stroud-based stochastic collocation methods [35] that may also outperform tensor-based solutions, the Small-Perturbation Automatic-Differentiation approach [36] that requires only small changes to existing deterministic codes, and uncertainty evaluation approaches based on the unscented transform [37] which are proven to be more efficient than MC techniques.The variety of existing solutions extends the range of stochastic problems that can be dealt with in an efficient framework and provides different lines of work that can be potentially combined with the efficiency of the FDTD computational scheme.