Conformally Mapped Polynomial Chaos Expansions for Maxwell's Source Problem with Random Input Data

Generalized Polynomial Chaos (gPC) expansions are well established for forward uncertainty propagation in many application areas. Although the associated computational effort may be reduced in comparison to Monte Carlo techniques, for instance, further convergence acceleration may be important to tackle problems with high parametric sensitivities. In this work, we propose the use of conformal maps to construct a transformed gPC basis, in order to enhance the convergence order. The proposed basis still features orthogonality properties and hence, facilitates the computation of many statistical properties such as sensitivities and moments. The corresponding surrogate models are computed by pseudo-spectral projection using mapped quadrature rules, which leads to an improved cost accuracy ratio. We apply the methodology to Maxwell's source problem with random input data. In particular, numerical results for a parametric finite element model of an optical grating coupler are given.

x y z dx dy Sketch of considered unit cell corresponding to the computational domain D. The blue arrow illustrates the incident wavevector k inc .

FIGURE 1
Scattering of periodic structure excited by an incident plane wave.
reported recently, see [12,13]. Although, these works equally rely on mapped Polynomial Chaos approximations, the transformations are linear (affine) and not based on conformal mappings. Also, the emphasis there is on high dimensional approximation instead of convergence acceleration.
The proposed numerical scheme is applied to quantify uncertainties via surrogate models for Maxwell's equations in the frequency domain. In particular, we consider the source problem on periodic domains with a plane wave excitation and uncertainties in the material interface geometry.
Such model equations can describe, for instance, the coupling into metal-insulator-metal (MIM) plasmon modes with subwavelength diffraction gratings, which is illustrated in Fig. 1a. Although illustrated by means of this particular application example, we note that the employed UQ methodologies apply in a much broader context. This paper is structured as follows: Section 2 contains a brief description of Maxwell's source problem. The uncertainty quantification part can be found in Section 3, where we briefly recall standard gPC before discussing the proposed extension based on conformal mappings. Section 4 reports numerical results for an analytical RLC circuit and the aforementioned optical grating coupler, before conclusions are drawn.

MAXWELL'S SOURCE PROBLEM
We consider Maxwell's source problem for periodic structures excited by an incident plane wave. For further details on this subject, we refer to [11,14]. We start with the time-harmonic curl-curl equation for the electric field phasor E in the computational domain D, where ω denotes the angular frequency, ε the complex permittivity and µr, µ 0 denote the relative and vacuum permeability, respectively. Note that (1) assumes absence of charges and source currents in D. Based on Floquet's theorem [14,Chapter 13], the computational domain D can be reduced to a unit cell of the periodic structure, as we assume a periodic excitation. Such a unit cell is depicted in Fig. 1b. Due to the oblique angle of the incident wave, the excitation has a different periodicity than the geometry and, hence, periodic phase-shift boundary conditions need to be imposed on the respective boundaries. To truncate the structure in the non-periodic direction, a Floquet absorbing boundary condition and a perfect electric conductor (PEC) boundary condition are applied. This leads to the boundary value where we refer to [11, Appendix A] for a derivation and definition of the functionals F (·), G(·).
We assume in the following that the complex permittivity ε depends smoothly on a parameter vector y ∈ Ξ ⊂ R N . These parameters can then be used to model variations in the refractive indices or extinction coefficients of the (different) materials in D, as well as changes in the geometry of the material interfaces inside the domain D. Following a standard Galerkin procedure, cf. [11], we then obtain a FE model in the form where ay(·, ·) is a continuous sesquilinear form, ly(·) is a continuous (anti)linear form and V denotes a discrete subspace of H (curl; D) [15], enforcing periodic phase-shift conditions on the traces at the periodic boundaries and homogeneous Dirichlet conditions at Γ z − . To achieve a curlconforming discretization of (7), we employ Nédélec's elements of the first kind [16] and 2nd order on a tetrahedral mesh of D. As QoI we consider the fundamental reflection coefficient Q(e(y)), i.e. a scattering parameter, which can be computed as an affine-linear functional of the electric field e in post-processing [11]. For brevity, we replace Q(e(y)) by Q(y) in the following.

UNCERTAINTY QUANTIFICATION
To account for uncertainty, we model the input parameters y as independent random variables (RVs) with joint probability density function ρ and image set Ξ ⊂ R N , where we assume in this section for brevity of notation that Ξ is given as the hypercube [−1, 1] N . Note that different image sets Ξ or stochastic dependence could also be considered, e.g. by a Rosenblatt transformation [17]. Additionally, we assume that the map Q : Ξ → C is holomorphic. Note that this assumption can often be justified for boundary value problems with random influences, see, e.g., [18]. Holomorphy of the solution of Maxwell's source problem with respect to general shape parametrizations was established in [19].
As discussed in the following, in this work we propose a method for surrogate modeling, where the basis functions are mapped polynomials based on gPC [7] combined with a conformal mapping. To compute the corresponding coefficients we rely on pseudo-spectral projection based on mapped quadrature rules [8].

Generalized Polynomial Chaos
For convenience of the reader, we briefly recall the standard polynomial chaos expansions, going back to Wiener [20]. Considering Gaussian random variables, any Q(y) with bounded variance, can be accurately represented using Hermite polynomials as basis functions. Employing the Askeyscheme [7], for different probability distributions ρ, basis functions Ψm : Ξ → R which are orthonormal w.r.t. the probability density ρ, i.e., can be obtained. We note that gPC can also be constructed for arbitrary densities ρ [21]. The gPC approximation is then given as where the sm ∈ C denote the gPC coefficients. In practice, in order to obtain a computable expression, the sum in (9) has to be truncated to M < ∞ and limited polynomial degrees are considered. The coefficients sm can then be determined in various ways, e.g. by regression or stochastic collocation, see [22] for an overview. Here we consider projection, i.e., The integral in (10) is usually not readily computable and is hence often approximated by numerical quadrature. Due to orthogonality of the basis, stochastic moments as well as variance-based sensitivity indices can then be calculated directly from the coefficients sm without further approximations, see [22]. These methods show spectral convergence, e.g., in the norm ||u|| L 2 ρ := E[u 2 ] [7]. In particular, if the map y → Q(y) is analytic, exponential convergence can be expected, as discussed in the following. Note that, for simplicity, we first consider the univariate case, i.e., N = 1, while generalizations to the multivariate case N > 1 will be discussed later.
We assume that Q 1D : [−1, 1] → C can be analytically extended onto an open Bernstein ellipse Er ⊂ C. A Bernstein ellipse Er is an ellipse with foci at ±1 and the size r is given by the sum of the length of semi-major and semi-minor axis. This is illustrated in Fig. 2a. Following [9], the error of the polynomial best approximation Q PC * M with degree M can be estimated as where · ∞ denotes the supremum-norm on [−1, 1] and the constant C B > 0 depends on the uniform bound of Q 1D in Er. Note that convergence in the supremum-norm implies convergence in the || · || L 2 ρ norm as well, as We further note that the additional aliasing error introduced by the discrete projection does not harm the convergence order for well-resolved smooth function, cf. [7, Chapter 3.6].

FIGURE 2
Illustration of conformal mapping approach.

Conformally Mapped Generalized Polynomial Chaos
Equation (11) shows that the convergence is connected to the region of analyticity, in particular the convergence order r depends on the size of the largest Bernstein ellipse not containing any poles of the continuation of Q 1D (in the complex plane). However, established procedures [23] inferring the regularity of parametric problems based on a sensitivity analysis, do not lead to elliptical regions, but rather prove analyticity in an -neighborhood of the unit interval. In this case, a conformal map g can be employed, which maps Bernstein ellipses to straighter regions and thus, enlarges the domain of analyticity, as illustrated in Fig. 2b. To this end, there are various mappings which could be employed, cf. [24]. Here, we focus for simplicity on the so-called 9-th order sausage mapping g(s) = 1 53089 (40320s + 6720s 3 + 3024s 5 + 1800s 7 + 1225s 9 ) introduced in [8], which represents a normalized Taylor approximation of the inverse sine function. Note that g maps the unit interval to itself, i.e., Conformal maps were employed in [8] to derive new numerical quadrature formulas, and have also recently been considered in the context of stochastic collocation methods [10,11]. In this work, we address the combination of conformal maps and polynomial chaos expansions. Based on the assumption that h := Q 1D • g has a larger Bernstein ellipse than Q 1D , and is hence better suited to be approximated with polynomials, we propose a new orthogonal basis whereΨm are orthonormal polynomials w.r.t. the transformed densitỹ We emphasize that {Φm} M m=0 forms an orthonormal basis w.r.t. the input probability distribution ρ. This can be shown by a change of variables y = g(s) where the last line holds by construction of the polynomialsΨm. Due to the orthogonality, the corresponding coefficients sm of the mapped can then be determined by the projection Note that, by abuse of notation, we use the same symbol sm for the gPC coefficients and the mapped gPC coefficients. The mapped polynomial best approximation Q * M converges as where h PC * M denotes the polynomial best approximation of h andr the size of a Bernstein ellipse Er on which an analytic continuation of h exists.
In particular, the convergence order of the mapped approximation Q M is given by the size of the largest Bernstein ellipse Er max which is fully mapped into the region of analyticity of Q 1D (y). Note thatrmax > rmax for any positive < 0.75 [11], and hence, a convergence improvement is to be expected in those cases, i.e., for functions analytic in such -neighborhoods. It should be mentioned nevertheless that this procedure does not always yield improved convergence rates. One can easily imagine poles located such that a Bernstein ellipse may lead to a larger region of analyticity than a strip-like geometry.
In the examples considered in this work, however, convergence acceleration could indeed be obtained.
To numerically compute (20), we derive mapped quadrature rules, cf. [8,24]. As pointed out in [9] for instance, Gaussian quadrature is derived from polynomial approximations and, hence, the convergence order also depends on the size of the Bernstein ellipse corresponding to the regularity of the integrand, see e.g. [8,Theorem 1]. Therefore, relying again the assumption that Q 1D • g has a larger Bernstein ellipse, we apply a change of ds.
The mapped quadrature scheme is then obtained by application of Gaussian quadrature w.r.t. the transformed densityρ 1D , i.e. quadrature nodes , to the transformed integrand in (22) sm ≈ Note that the mapped quadrature nodes are obtained asŷ (i) := g(ỹ (i) ), while the mapped weights are given asŵ (i) :=w (i) . Due to (14), it is ensured that the mapped quadrature nodesŷ (i) do not require the evaluation of the analytic continuation of Q 1D in the complex plane, which is, in practice, not always possible. A convergence improvement is expected based on the assumption that the transformed integrand in (22) has a larger Bernstein ellipse. For further details on mapped quadrature schemes, we refer to [8]. However, we note the (minor) differences that in this work we employ Gaussian quadrature w.r.t. the transformed densityρ 1D to derive the mapped quadrature scheme, while [8] only considers unweighted Gaussian quadrature and, thereby, takes g (s) as part of the integrand (instead of the weight).
We proceed with a discussion of the multivariate case N > 1. To this end, we introduce the multivariate mapping g(s) = [g 1 (s 1 ), . . . , g N (s N )].
In this work, we employ, for simplicity, the same mapping (13) for all parameters, i.e. g 1 = . . . = g N = g. However, different choices would be possible as well. We also note that, for the trivial mapping g triv : s → s standard polynomial chaos expansions would be recovered. For each parameter y i with univariate probability density function (PDF) ρ i , we define the transformed PDFρ i (y i ) := ρ i (g i (y i ))g i (y i ). The corresponding transformed joint PDF is then given byρ(y) =ρ 1 (y 1 ) . . .ρ N (y N ). In the following, we denote by {Ψm}m an orthonormal polynomial basis w.r.t. to the transformed densityρ, i.e.
The coefficients of the multivariate mapped approximation where we consider for simplicity a tensor-product construction of maximum degree p, can then again be obtained by projection To evaluate the multi-dimensional integral in (27), we employ mapped Gaussian quadrature. In this case the mapped nodes and weights are given bŷ y (i) := g(ỹ (i) ) andŵ (i) :=w (i) , respectively, where, in turn,ỹ (i) andw (i) are the nodes and weights of a Gaussian quadrature w.r.t.ρ.
Finally, we emphasize that, since the mapped representation (26) uses an orthogonal basis, the coefficients sm can be used to directly compute stochastic moments as well as variance-based sensitivity indices. For instance, the mean value is given by where we employed, that the mapped basis function Φ 0 is constant on Ξ, as well as the orthonormality condition (24). Accordingly the variance is given by  Additionally, Sobol sensitivity indices [25], based on a decomposition of the variance, can also be directly derived from the coefficients. Regarding the estimation of Sobol indices, we will focus on the so-called main-effect (1st order) and total-effect (total order) indices. We define the multi-index Then, the main-effect and total-effect Sobol indices, S main n and S total n , respectively, are given as

APPLICATION
We apply the UQ methods presented in the last section to two model problems. We first consider an academic example of an stochastic RLC circuit, since there is a closed-form solution available which allows us to illustrate the main ideas of the proposed approach in detail. We then consider the optical grating coupler [2], which is a non-trivial benchmark example from nanoplasmonics.

RLC circuit
We consider the model of an RLC circuit, as illustrated in Fig. 3a. Assuming harmonic time dependency, the electric current i is given by We consider, arbitrarily chosen, an angular frequency ω = 10 4 s −1 , exciting voltage ue = 1 V, capacitance C = 10 µF, and a (rather small) resistance of R = 1 Ω. Additionally, we consider a variable inductance L(y) = 1 mH + 0.25 mH · y. The parameter y is then modeled as a uniformly distributed random variable with probability density function ρ = U (−1, 1), such that a stochastic model is obtained. As QoI Q, we consider the amplitude of the current Q := |i|. Fig. 3b shows the parametric dependency of the QoI |i| with respect to y, which is analytic for y ∈ [−1, 1].
However, the continuation in the complex plane has poles at This complex conjugate pole pair limits the size of the largest Bernstein ellipse, where Q(y) is analytic, which is illustrated in Fig. 4a for different values of R.
In each case, we compute gPC approximations of increasing order for Q(y) using the Chaospy toolbox [26]. In particular, the gPC coefficients of an M −th order approximations are computed by pseudo-spectral projection using Gaussian quadrature of order M + 1. The accuracy of the surrogate models is then quantified in the empirical L 2 ρ norm. In particular, we apply cross-validation using N cv = 1000 random parameter realizations y (i) cv drawn according to the probability density ρ, to compute the error        Fig. 4b and Fig. 4c, respectively. The plots confirm (11) numerically, showing a decreasing convergence order for decreasing values of R corresponding to decreasing sizes of the associated Bernstein ellipses. Note that, according to (35) a similar behaviour as for decreasing damping can be expected for increasing amplitudes of the considered input variation.
Next, we apply the conformally mapped gPC expansions proposed in the last section. The implementation is done in Python based on Chaospy [26]. Fig. 5a shows the transformed density (16) for a uniform input distribution ρ. Fig. 5b depicts some exemplary basis functions of gPC and mapped gPC. Note that the gPC basis functions are in this case Legendre polynomials, while the mapped basis functions, given by (15), are no polynomials.
We then study the convergence of the corresponding surrogate models, where mapped quadrature of order M + 1 is used to compute the mapped gPC expansions of order M . Fig. 5c, Fig. 5d and Fig. 5e demonstrate the improved convergence order of the mapped approach, in terms of the cross-validation error, as well as the accuracy of the computed mean value and the computed standard deviation.

Optical Grating Coupler
We now consider the FE model of an optical grating coupler [2], which was introduced in the beginning, see Fig. 1a. The structure's design [27] is shown in Fig. 6. A plane wave at optical frequency hits the surface of the grating coupler. The incident wave couples with a MIM plasmon mode, which propagates along the metallic surface. It is found that the MIM resonance has a significant shift (in energy) as a function of the grating depth [2] and therefore, it is of great interest to evaluate the influence of nano-technological manufacturing imperfections.
We use FENICS [28] for the discretization and implement a design element approach [29] for the geometry parametrization. The numerical model is described in greater detail in [11]. Note that we only consider periodic variations, modeling a systemic offset in the fabrication process, and do not address local uncertainties leading to different unit cells. Readers interested in the latter case are referred to [30]. The fundamental scattering parameter is considered as QoI Q ∈ C. We consider three sensitive geometrical parameters as uncertain, in particular the thicknesses of the upper gold layer t 1 = 12 nm + ∆y 1 , the thickness of the dielectric layer t 2 = 14 nm + ∆y 2 and the grating depth T = 20 nm + ∆y 3 , as illustrated in Fig. 6. We model those parameters as independent beta distributed RVs in the range of ±∆ = ±2 nm. The corresponding shape parameters are chosen such that a normal approximation is approximated. The corresponding probability distribution ρ i of the RVs y i , i = 1, . . . , 3 is shown in Fig. 7a, together with the transformed densityρ i . The univariate gPC polynomials which are Jacobi polynomials in this case, as well as the mapped polynomials are illustrated in Fig. 7b.

Decay of Fourier Coefficients
We first study the decay of polynomial coefficients to numerically investigate the smoothness of the mapping from the input parameters to the complex S-parameter Q and justify the use of (mapped) polynomial approximations. It has been shown, see e.g. [31,Lemma 2] where Legendre polynomials are considered, that if this mapping is smooth, the Fourier coefficients sm of an N −variate gPC approximation decay exponentially, i.e.
It can be seen that the maximum Fourier coefficient is expected to decay exponentially with an increasing maximum-degree w.
We construct a multivariate gPC approximation with a tensor-product basis of order mmax = 15. The multivariate integrals of the pseudospectral projection are then computed by a Gauss quadrature of order 17. All coefficients sm are plotted in Fig. 8 in red color, where an exponential decay can indeed be observed. This can be seen as a numerical indicator for smoothness of the approximated mapping Q(y). Additionally, we also construct a mapped approximation of same order and plot the corresponding coefficients in black color. It can be observed that the mapped coefficients exhibit a faster convergence and hence the mapped approach can be expected to show, again, an improved convergence.

Uncertainty Quantification
Next, we consider approximations of the magnitude of the S-parameter |Q(y)| using (mapped) tensor-product gPC expansions of increasing order M , where pseudo-spectral projections of order M + 1 is employed to compute the coefficients. Fig. 9a compares gPC and the proposed mapped counterpart in terms of the L 2 ρ -error (36), in particular, again, by cross-validation with 10 3 random parameter realizations. It can be observed that the mapped approach converges about 30% faster w.r.t. the order M than gPC. However, the respective computational gain grows, in this case, exponentially w.r.t. the number of inputs and, hence, the required number of model evaluation to reach a prescribed accuracy reduces roughly by a factor of 2. Similar findings hold for the stochastic moments, in particular, we present the convergence of the mean value in Fig. 9b and the computed standard deviation in Fig. 9c. In this case, the reference solutions are obtained by Gaussian quadrature of order 30.  Fig. 10. The thickness of the dielectric layer t 2 is identified as the most influential parameter. We note that there is a significant difference between the main-and total-effect indices. In particular, the sum of the first order indices is only 34%, while the remaining 66% can be attributed to strong coupling effects among the parameters.

CONCLUSIONS
In this paper an efficient surrogate modeling technique for quantifying uncertainties in the material and geometry of high-frequency and optical devices was presented. The proposed method is based on gPC to achieve spectral convergence. Through a combination with conformal maps we were able to enlarge the region of analyticity. This lead to an improved convergence rate, which was numerically demonstrated for two benchmark problems. In particular, the approach showed significant gains in either accuracy or computational cost, without adding any relevant extra computational effort. Due to orthogonality of the proposed basis, stochastic moments as well as Sobol indices can be directly obtained from the coefficients.
It is worth noting that this technique can also be combined with other techniques for convergence acceleration such as adjoint-error correction, sparse-grids and (adjoint-based) adaptivity for the multivariate case [11].