Uncertainty quantiﬁcation of silicon photonic devices with correlated and non-Gaussian random parameters

: Process variations can signiﬁcantly degrade device performance and chip yield in silicon photonics. In order to reduce the design and production costs, it is highly desirable to predict the statistical behavior of a device before the ﬁnal fabrication. Monte Carlo is the mainstream computational technique used to estimate the uncertainties caused by process variations. However, it is very often too expensive due to its slow convergence rate. Recently, stochastic spectral methods based on polynomial chaos expansions have emerged as a promising alternative, and they have shown signiﬁcant speedup over Monte Carlo in many engineering problems. The existing literature mostly assumes that the random parameters are mutually independent. However, in practical applications such assumption may not be necessarily accurate. In this paper, we develop an efﬁcient numerical technique based on stochastic collocation to simulate silicon photonics with correlated and non-Gaussian random parameters. The effectiveness of our proposed technique is demonstrated by the simulation results of a silicon-on-insulator based directional coupler example. Since the mathematic formulation in this paper is very generic, our proposed algorithm can be applied to a large class of photonic design cases as well as to many other engineering problems.


Introduction
The silicon photonic technology has emerged as a promising alternative to electrical interconnects due to its ability to achieve higher bandwidth and lower power dissipation [1]. Thanks to its ease of integration with the CMOS process, silicon photonics can remarkably reduce fabrication costs, increase the integration scale and improve the overall system performance [2]. However, silicon-based optical devices are very sensitive to fabrication process variations, leading to potentially significant device performance degradations and potential system failures [3][4][5][6][7][8].
There is a lack of uncertainty quantification techniques for silicon photonics, however some results have been reported recently for nanometer integrated circuits (IC) [9][10][11][12] and microelectromechanical systems (MEMS) [13,14]. The mainstream stochastic simulation technique in commercial design software is Monte Carlo [15]. Despite some recent improvements [16], Monte Carlo simulation is still inefficient for most design problems. This limitation is clearly evident in silicon photonics where reliable Monte Carlo results can be obtained only by performing a large number of simulations, resulting in a very high computational cost.
Generalized polynomial chaos (gPC) expansion techniques have been proposed [17] to efficiently approximate a stochastic solution dependent on a set of random parameters. Furthermore, they have recently been exploited in the analysis of photonic circuits [18]. Based on such expansions, fast stochastic spectral methods have been developed [19][20][21]. Standard spectral methods include an intrusive method we refer to as stochastic Galerkin (SG) [19], a nonintrusive method called stochastic collocation (SC) [20,21], and their hybrid variant called stochastic testing [9][10][11]. In most existing publications [18,[22][23][24][25], the input parameters are assumed to be mutually independent, or Gaussian-correlated, in which case one can convert the random variables to a set of independent ones by a linear transform. However, in practical applications such assumption may not be necessarily accurate.
In this paper, we propose a computational framework to quantify the uncertainties in silicon photonic devices caused by non-Gaussian correlated random parameters. Our algorithm is based on the theoretical results of [26] which has been rarely utilized in practice due to relevant issues faced in its numerical implementation. In our work, we assume that the process variations are described by a Gaussian Mixture (GM). Gaussian mixture is quite a reasonable modeling technique for process variations, since it can capture both correlated and uncorrelated process variations. Based on this assumption, we first develop a numerical scheme to construct a set of orthogonal basis functions to represent the stochastic solution. We further propose an implementation of stochastic collocation to compute each of the coefficient associated with each of the basis functions. Differently from our preliminary conference paper [27], this paper provides also all the theoretical and numerical details for computing the basis functions and for implementing stochastic collocation with correlated random parameters. As an example, a silicon-on-insulator (SOI) based directional coupler is considered and the impact of tolerances on geometric dimensions is evaluated. Moreover, the accuracy of our proposed algorithm is compared with Monte Carlo in our numerical experiments. This paper is structured as follows. In Section 2, we give a brief review of orthogonal basis functions with respect to a joint probability density measure, as well as some basics about Gaussian mixture models. Section 3 consists of two parts describing our proposed computational techniques. The first part demonstrates how to construct the orthogonal bases given a Gaussian-mixture density function; the second part presents how to implement efficiently a stochastic collocation scheme to obtain the stochastic solution. In Section 4, we report the simulation results of a SOI based directional coupler under process variations. Finally, Section 5 concludes this paper.

Background review
This section briefly reviews some background related to orthogonal basis functions [26] and Gaussian mixture models.

Orthogonal basis function construction
Let ⃗ ξ = {ξ 1 , · · · , ξ N } ∈ R N be N random variables with a given joint probability density function p ⃗ ξ on the support Ξ: We further denote the marginal probability density function of the i-th random variable by p i (ξ i ): Assume that {φ (i) k (ξ i )} are a set of generalized polynomial-chaos bases for ξ i [17], where k denotes the degree of a polynomial. Such basis functions are orthogonal with respect to the marginal density function p i : where ⟨⟩ p i denotes an inner product defined by the marginal probability density function p i (ξ i ).

Gaussian mixture
In this paper, we assume that a Gaussian-mixture model has already been constructed based on some measurement data, using for instance the Expectation-Maximization (EM) algorithm [28]. We focus here instead on how to predict the effect on the performance of a silicon photonic device influenced by process variations. In practice, there are correlated and uncorrelated process variations. In this paper, we aim at developing numerical techniques that can deal with correlated process variations, since there are already many techniques available for simulating design problems with uncorrelated process variations. A Gaussian mixture with M mixed terms is the weighted sum of M multivariate Gaussian densities. Purely for simplicity of exposition of mathematical derivations, we set here N = M = 2. However, all results presented in this paper are valid and can be extended in a straightforward way to the more general cases where N and M are not equal to 2 and N is not necessarily equal to M.
With a Gaussian mixture model, the process variations can be described by the following distribution with where ⃗ ξ can be the geometrical or material parameters of a silicon photonic device. It is wellknown that the marginal probability density function of a multivariate Gaussian is a univariate Gaussian. Therefore, the i-th random parameter follows a univariate Gaussian-mixture distribu- and its marginal probability density function is

Uncertainty quantification with non-Gaussian correlated parameters
In this section, we first develop some numerical techniques to build the orthogonal basis functions introduced in Section 2. We then describe how to implement efficiently a stochastic collocation technique to calculate a stochastic quantity of interest that depends on some correlated non-Gaussian process variation parameters. This section provides all mathematical details to construct a stochastic solution based on stochastic collocation, and the flows are summarized in Algorithm 1 and 2. For readers who are more interested in the applications of the proposed computational technique, please refer to Section 4 for a directional coupler example.

Basis function computation
Let u(x, ⃗ ξ ) be the quantity of interest, smoothly dependent on the process variation parameters.
Here u can be for instance the effective phase index, the resonant wavelength or the power dissipation of a silicon photonic device. Vector ⃗ ξ could contain for instance the dimensions and/or the material properties of a device, and x can be the frequency or wavelength. For example, in order to investigate the uncertainties caused by geometric variations of a directional coupler, we can set u to be the power coupling coefficient, and ⃗ ξ to be the exterior and interior waveguide sidewall position. Meanwhile, x can be set as the operating frequency or wavelength that also influences the power coupling coefficient. Given the basis functions Ψ ⃗ α ( ⃗ ξ ) described in Section 2, we can approximate u as when u has a bounded variance.
, ⃗ α is a multiindex which consists of N-tuples, and a full-tensor truncation is used for gPCe. The basis functions are still not known, and thus our first task is to compute them given the Gaussian-mixture joint density function. In order to construct Ψ ⃗ α , we first need to compute the univariate basis functions φ (3) for each random parameter ξ i . The univariate basis functions can be built using the three term recurrence relations [29], and the construction of such univariate orthogonal polynomials is shown in the following. Since our numerical scheme is applicable to all random parameters ξ i , the parameter index i in φ (i) k (ξ i ) and p i (ξ i ) will be omitted in the following derivation.
The three term recurrence relations is In order to ensure that {φ k } k≥0 is an orthogonal set with respect to measure p, the coefficients a k and b k are defined as In order to compute a k and b k , we need to evaluate with p(ξ ) defined in Eq. (7). It is non-trivial to obtain Eq. (12) in an analytical way. However, , the integral in Eq. (12) can be rewritten as As a result, a Gauss-Hermite quadrature rule can be used to numerically compute Eq. (13) with high accuracy. For instance, the two terms in Eq. (13) can be computed as where x j and w j are the j-th well-known Gauss-Hermite abscissa and weight respectively [30] , and q is the number of quadrature points. The inner products ⟨φ k , φ k ⟩ p and ⟨φ k−1 , φ k−1 ⟩ p can be calculated in a similar way: With a k and b k computed with the above procedures, the formula in Eq. (9) can be used to calculate the basis function for each parameter ξ i . Finally, the multivariate basis function Ψ ⃗ α can be obtained according to Eq. (4).

Implementation of stochastic collocation
The coefficients C ⃗ α (x) in Eq. (8) can be computed based on stochastic collocation, and the implementation is detailed below.

Benchmark description
Due to the high refractive index contrast between silicon and silica, silicon photonic devices are extremely sensitive to variations in the geometry, for instance in waveguide width. Such geometrical variations can cause significant fluctuations in the effective phase index (n eff ) and can further lead to performance degradations in many photonic devices such as directional couplers and ring resonators. In practical fabrications, due to tolerances of the fabrication process, the fabricated waveguide width (denoted as W ) and the gap (denoted as g) between two waveguides of a directional coupler (as shown in Fig. 1) are different from chip to chip and wafer to wafer. Here, we assume that exterior and interior waveguide sidewall position, ∆W e and ∆W i , are correlated random parameters. The two coupled waveguides are assumed identical for simplicity and ∆W i is assumed to be different from ∆W e because of the typical proximity effects in the etching process. The fabricated W and g are related to the originally designed width (W 0 ) and (g 0 ) as output phase difference. The power coupling coefficient K(z) is a sinusoidal function of the coupler's length z where δ is the field coupling coefficient and it is equal to the effective phase index difference between the symmetric mode and the asymmetric mode of a coupler times a factor π/λ , where λ is the wavelength of concern. On the other hand, the output phase difference is always π/2 due to the assumed symmetry. Note that exact density functions for ∆W e and ∆W i can be extracted by measuring a large set of practical chips. However, in this paper, for the purpose of algorithm verification, we assume ∆W e and ∆W i are correlated, since there are already extensive techniques for the non-correlated process variations. Under the above assumption, our computational methods developed in Section 3 can be employed to analyze how such geometric uncertainties affect the power coupling coefficient. According to Eq. (24), it will be more convenient to approximate δ using gPC, since we can compute K(z) for any length z with δ at hand. Specifically, we approximate δ (λ , ∆W e , ∆W i ) by a linear combination of a set of orthogonal functions {Ψ ⃗ α (∆W e , ∆W i )}, where λ is the wavelength. When computing C ⃗ α (λ ), we need to evaluate δ for each multidimensional quadrature point in Eq. (20). For each sample, a finite difference mode solver is called to solve for δ with the given value (∆W e , ∆W i , λ ) using Eq. (23). In our numerical experiments, we focus on synchronous couplers where the nominal width W 0 is 0.4 µm and the nominal gap g 0 is 0.2 µm. The waveguides height are fixed to 0.22 µm, the refractive index of silicon and silica are 3.48 and 1.445 respectively, and the wavelength of interest is 1.55 µm. The simulated δ with nominal width and gap is 0.11. All simulations are performed using a personal computer with a intel-xeon CPU X5650 and 24G of RAM.

Numerical results
In order to validate our computational methods, we use the following Gaussian mixture for simulation: With a proper Gaussian mixture, we can approximate any measured distribution. For simplicity we have included only two terms in the Gaussian mixture here. The numbers in the mean and covariance matrix are picked to have the total variations of coupler's width (∆W ) and gap (∆g) range from around 0 to 25 nm, as shown in Fig. 2. The chosen numbers are reasonable numbers in current fabrication process tolerance. To show our algorithm's capability of dealing with correlated parameters, the weights a and b are chosen close to 0.5, and the exact values are 0.6 and 0.4 respectively. Note also that the choice of different distributions for ∆W e and ∆W i here allows taking into account the proximity effect that occurs in the gap during the fabrication.
Our algorithm constructs an approximation of the parameter-dependent δ in the form of Eq. (8), and then computes K(z) using Eq. (24). By doing so, we can easily evaluate its density function without calling the expensive mode solver every time in Monte-Carlo simulation. Figure 3 shows the probability density function of δ at λ = 1.55 µm. The solid line represents δ obtained from our stochastic collocation with level ⃗ m = (5,5). Here ⃗ m = (5, 5) means that a 5-level Gauss-Hermite quadrature rule is used for both parameter dimensions to evaluate Eq. (20). For parameter ξ i , the level m i is determined as where q i is the number of quadrature points used for ξ i . Therefore, 9 quadrature points are used in each dimension to obtain the results in Fig. 3, i.e. q 1 = q 2 = 9 in Eq. (20). Note that to avoid aliasing error in the gPC coefficient, the relation between ⃗ m and ⃗ t in Eq.
In order to verify the accuracy of our approach, the simulated probability density function from standard Monte Carlo with 10 4 samples is also plotted as a dashed line in Fig. 3. The result from our stochastic collocation is consistent with that from Monte Carlo method. The nominal value of δ for W 0 = 400 nm and g 0 = 200 nm is 0.11, while the analysis shows an expected value asymmetrically distributed around δ = 0.102. With the pdf of δ at hand, we can compute K(z) using Eq. (24), and the result is shown in Fig. 4. Figures 4(a) and 4(b) show the simulated probability density functions of the power coupling coefficient at z = 7.1376µm and z = 2.9250µm, which are the lengths that correspond to K 0 = 0.5 and K 0 = 0.1 respectively. K 0 is the nominal power coupling coefficient for the nominal width W 0 = 400 nm and gap g 0 = 200 nm. From Fig. 4, it is clearly seen that K deviates from its original desired values under process variations. Notice that with the choice of the parameters in our example (W 0 = 400 nm and g 0 = 200 nm), the coupler coefficients are on average smaller than designed (K 0 = 0.5, and K 0 = 0.1). This can be explained by looking at the average waveguide width and gap in our example: the waveguide width are always larger than the nominal values (around 15 nm in average) while the gap is always smaller (around 12.5 nm in average). For example, the computed K(z = 7.1376µ m) with W = 415 nm and g = 187.5 nm is around 0.44, which is smaller than the designed value 0.5 and complies the simulated results in Fig. 4(a). Similar results can be obtained for the case of K(z = 2.9250µ m) with W = 415 nm and g = 187.5 nm.
In order to show the convergence of our algorithm, we gradually increase the level numbers and compare the truncation error of δ (λ , ∆W e , ∆W i ). We use the 10 4 samples Monte Carlo simulation result as the reference solution and plot the L 2 -norm of the error as a function of the level parameter. The result in Fig. 5 shows that the error decreases when more basis functions    and more Gauss quadrature points are used. In addition, to quantitatively show the accuracy of our solver, Table 1 compares the mean value and standard deviation of the simulated δ from Monte Carlo and from our technique. It clearly shows that Monte Carlo result converges towards that from our proposed solver as the number of samples increases. Moreover, to achieve a similar level of accuracy in the mean and standard deviation of δ , Monte Carlo uses 10 4 samples and 97 minutes of CPU time, whereas our method consumes about 1 minute and 45 seconds when 9 quadrature points are used in each dimension. This clearly shows the superior efficiency of our technique over Monte Carlo simulation, and the speedup factor is about 55. Remark: In our implementations, we have observed slower convergence rates for some cases even when the δ is smoothly dependent on the process variations parameters ∆W e and ∆W i . The slow convergence rate is caused by the theoretical limitations of the basis functions proposed in [26]. Observing Eq. (4), we note that the basis functions Ψ ⃗ α ( ⃗ ξ ) might be highly non-linear when the parameters are correlated. Thus, for some cases many basis functions are needed to approximate a smooth quantity of interest (which is δ in our example), resulting in a slower convergence rate for the solution. In order to improve the convergence rate, we are working on a whole new set of basis functions which can potentially achieve high accuracy with much fewer basis functions.

Conclusions
In this paper, we have developed a stochastic collocation scheme to simulate silicon photonic devices with non-Gaussian correlated parameters described by a Gaussian-mixture joint probability density function. The numerical techniques for constructing a set of orthogonal basis functions and for implementing stochastic collocation are presented. The proposed numerical