We demonstrate the use of stochastic collocation to assess the performance of photonic devices under the effect of uncertainty. This approach combines high accuracy and efficiency in analyzing device variability with the ease of implementation of sampling-based methods. Its flexibility makes it suitable to be applied to a large range of photonic devices. We compare the stochastic collocation method with a Monte Carlo technique on a numerical analysis of the variability in silicon directional couplers.
© 2016 Chinese Laser Press
Integrated photonics, and in particular, silicon photonics, is rapidly enabling complex photonic functions on a chip . However, variability due to fabrication processes and operational conditions limits the complexity of the circuits that can be implemented [2,3]. A proper analysis of the effects of the variability in geometrical, electrical, and optical parameters on the performance of the photonic building blocks and circuits has become crucial. Indeed, the variability introduced by the variations of the manufacturing process is a primary source of degradation of larger photonic circuits, especially when wavelength-selective filters are implemented . The primary functional parameters that are affected are the waveguide propagation constants (the effective index) and the coupling coefficients in coupling structures. The Monte Carlo (MC) method  is considered the standard approach for variability analysis, thanks to its accuracy and ease of implementation. Unfortunately, the MC analysis has a slow convergence rate, and it requires a large number of data points (simulations or measurements). Therefore, MC has a high computational cost, considering that accurate simulations of photonic devices can be time and resource intensive.
The generalized polynomial chaos (gPC) expansion has been applied in several domains as an efficient alternative to the classic MC method [6–8], and, recently, it has been proposed for the variability analysis of photonic devices [9,10].
The gPC-based modeling approach aims at expressing a stochastic process as a series of orthogonal basis functions with suitable coefficients and gives an analytical representation of the variability of the system on the random variables under consideration .
In this paper, we propose a stochastic collocation (SC) method as an efficient alternative to characterize photonic devices under the effect of uncertainty. The fundamental principle of the SC approach is to approximate the unknown stochastic solution by interpolation functions in the stochastic space. The interpolation is constructed by repeatedly solving (sampling) the deterministic problem at a predetermined set of nodes in the stochastic space. This approach offers similar high accuracy and efficiency as the stochastic gPC method, but, at the same time, it is easy to implement, such as sampling-based methods (e.g., MC approach).
We apply the SC method to analyze the variability of a key building block for silicon photonic circuits: the directional coupler (DC). This device is essential in the construction of wavelength filters, as it implements an arbitrary fractional power coupling, but, at the same time, it is extremely sensitive to fabrication variations: a small shift in linewidth or thickness of the core can dramatically change the coupling coefficients.
The paper is organized as follows: Section 2 gives a general introduction of SC methods. It provides the essential mathematics knowledge that readers need to know to understand its application in the photonic domain. Section 3 uses a DC as an example to test the performance of SC in performing variability analysis of photonic devices. Section 4 draws the conclusions.
2. STOCHASTIC COLLOCATION METHODS
SC methods are based on interpolation schemes to compute stochastic quantities. The interpolation is constructed by repeatedly solving (sampling) the deterministic problem at a predetermined set of nodes in the stochastic space  (also defined as collocation points). Indeed, a stochastic process can be expressed as
For a photonic device, the process could correspond to the functional parameters such as the waveguide propagating constants and the coupling coefficients in coupling devices. The stochastic variables correspond to device properties affected in a stochastic way by fabrication and operational conditions (e.g., waveguide line, width, or temperature).
In Eq. (1), different types of interpolation schemes can be adopted (e.g., piecewise linear [12,13], Lagrange [11,14], or interpolation methods that belong to the general class of positive interpolation operators such as multivariate simplicial methods ). However, the key issue for this approach is the selection of the support nodes, such that when using the minimal number of nodes one achieves good approximation.
For example, if the Lagrange interpolation scheme is chosen, the element in Eq. (1) for a 1D interpolation can be expressed as1) becomes 3) is 4), the number of nodes required by the full tensor product increases rapidly with the number of random parameters . For example, if three random variables are considered and 10 collocation points are used for each parameter, a total of 1000 nodes are required by the full tensor product approach. Hence, the performance of the photonic device under study must be evaluated for 1000 different combinations of the random variables considered, leading to an expensive computational time.
The required number of nodes can be significantly reduced by adopting sparse grids in the stochastic space, based on the Smolyak algorithm [12,16–19]. By choosing the collocation points correctly, the Smolyak algorithm drastically reduces the total number of nodes used in the interpolation with respect to the full tensor product approach while preserving a high level of accuracy.
It is important to remark that the SC models are expanded using interpolation functions of independent random variables . In the general case of correlated random variables, decorrelation can be obtained via a variable transformation, such as the Nataf transformation  or the Karhunen–Loéve expansion .
The stochastic moments (mean, variance, …) can be computed utilizing analytical formulas and then efficiently, once the analytical form of the interpolation functions has been decided. For example, if the random variables are defined in the sample space , the mean of is defined as1) in Eq. (5) leads to 6), an efficient numerical solution can be used (e.g., by MC analysis of the interpolation model or numerical integration). Finally, it is important to remark that it is not possible to define a priori the speed-up of a generic SC modeling technique compared with the MC method. Indeed, the number of nodes needed to compute an accurate SC model (which is directly related to the efficiency of SC methods, as described above) cannot be decided upfront because it depends on the following factors:
- • the impact of the chosen random variables on the variations of the stochastic process considered (dynamic stochastic processes require a higher number of collocation points);
- • the interpolation scheme adopted (the more powerful the interpolation scheme, the fewer nodes are needed);
- • the sampling strategy adopted (efficient sampling strategy limit the number of collocation points used);
- • the number of random variables considered (the higher the number of variables the more collocation points are needed).
However, it has been proven in the literature that, for a limited number of random variables (indicatively less than 10) SC methods are much more efficient with respect to the MC analysis (see [11,12,14]). For stochastic processes depending on a high number of random variables, the efficiency of SC methods is significantly reduced.
Two approaches can be used to increase the efficiency of an SC modeling technique. Using nested sampling schemes allows us to adaptively choose the collocation points (additional details are provided in Section 3.E and Appendix 5.B). Adopting adaptive sparse grids  reduces the nodes requirement, which is especially useful when a high number of random variables is considered. For a more detailed reference on SC methods, we refer the reader to [11,12,14,16].
3. DIRECTIONAL COUPLER EXAMPLE
A. Benchmark Description
We demonstrate the use of SC for integrated photonics through the analysis of a DC in a silicon photonics platform. Silicon photonics is rapidly gaining adoption due to its potential for large-scale integration and volume manufacturing. . However, the same high-index contrast that enables dense integration also makes silicon photonic waveguides extremely sensitive to small imperfections in their geometry. Also, the high thermo-optic coefficient of silicon makes silicon photonic devices temperature sensitive. For instance, a change in linewidth and thickness of waveguide would noticeably vary the effective index () of optical modes, resulting in a shift of 1 nm in the response of a wavelength-selective filter. This variation leads to performance degradation in devices such as DCs, Mach–Zehnder interferometers (MZIs) and ring resonators, and limits the number of devices that a single circuit can have.
Power coupling devices are essential parts for splitting and combining power in photonic circuits. Among them, DCs are widely used for their simplicity in layout and easy-to-understand operation. An advantage of DCs compared with other couplers is that the coupling ratio can be continuously adjusted by choosing the length of the coupling section. Furthermore, DCs constitute the building blocks of many larger photonic devices such as rings, MZIs, and so on.
A DC consists of two parallel waveguides and connecting bend sections: the light in a single waveguide is mostly confined in the silicon core, but an exponentially decaying field extends into the cladding. When two waveguides are brought in proximity, the modes of the two waveguides couple and form two supermodes with opposite symmetry (an even and an odd mode). The beating of these supermodes translates in a sinusoidal power transfer from one waveguide to the next and back along the propagation axis . The power coupling in a DC can be expressed as 
The rate of coupling is defined as the field coupling coefficient , which is determined by the geometry of the coupler cross section, such as the waveguide width, thickness, and gap between the waveguide cores.
Let us assume that, for simplicity, the two waveguides in a DC are identical. As a result, the straight section of the DC layout is defined by three parameters: the waveguide width , thickness , and gap (Fig. 1). Furthermore, we assume that, in the lithography process, the centers of the two waveguides are located at the designed position. It is a good assumption for optical lithography techniques but might be less accurate for e-beam written devices. With this assumption, the sum of the gap and the half-waveguide width is constant, as shown in Fig. 1. Therefore, in our example, we can describe the full geometry of the DC with only two parameters: and .
In this study, we will use the SC technique to find out how geometry variability influences the DC performance, namely, the coupling coefficient . Indeed, due to the fabrication variations, the fabricated linewidth , thickness , and gap are different on the value chosen during the design phase. To prove the robustness and modeling power of the proposed approach, we assume the width and thickness of the DC as correlated random variables, rather than independent, following the Gaussian distribution. It is not an unrealistic assumption: thickness variations could induce over-etching on the sidewalls.
It is good to note that the SC methods can deal with random variables with arbitrary distribution. It is therefore not necessary that the and adhere to a Gaussian distribution.
B. Simulation Setup
According to the theory of supermodes, we can write the coupling coefficient as 
To calculate of a given geometry, we define the DC structure accordingly and simulate and in the mode solver Fimmwave using its film matching mode (FMM) solver. For later performance comparisons, all simulations are performed on a computer with an Intel Core i5 2500 quad-core CPU clocked at 3.3 GHz and 8 GB of memory.
C. Problem Definition
As mentioned in Section 3.A, we considered the coupling coefficient of a DC as a stochastic process depending on two correlated random variables with Gaussian PDFs: the width and thickness . Hence, the joint PDF of the two random variables considered is defined as
To validate the robustness of the proposed method, and are chosen equal to 2% and the correlation coefficient , which is a challenging example to study because the coupling coefficient is quite dynamic with respect to the parameters considered (see Fig. 2). The proposed method is discussed in detail in the following and summarized in Fig. 3.
D. Variable Transformation
Now, SC methods in the form of Eq. (1) deal with independent random variables, as described in Section 2. Hence, to fit the problem into the SC framework, first of all, it is necessary to express the coupling coefficient in two independent Gaussian random variables, starting from the correlated random variables , defined by Eq. (9). As mentioned in Section 2, such decorrelation can be obtained via a variable transformation. Thanks to the Karhunen–Loéve expansion , it is possible to express the vector of correlated Gaussian random variables in the vector of uncorrelated Gaussian random variables with zero mean and unit variance as
E. SC Model Computation
To compute an SC in the form of Eq. (1), the first step is choosing the interpolation scheme: the Lagrange interpolation scheme is adopted in this example for its modeling power and ease of implementation. Next, a rule that guarantees a good quality of the approximation must be used to choose the collocation points for each random variable: and . In this example, the Clenshaw–Curtis rule is adopted : the collocation points for each random variable are the extrema of the Chebyshev polynomials. Now, the total number of nodes could be obtained using the full tensor products of the nodes chosen for each random variable, but it would not be efficient, as discussed in Section 2. Instead, the nodes are chosen over a sparse grid based on the Smolyak algorithm. Indeed, the adoption of the Smolyak algorithm allows building our SC model by using only a subset of all the collocation points given by the full tensor product . Furthermore, the collocation points chosen by the Smolyak algorithm based on the Clenshaw–Curtis rule are nested: if additional nodes are required to accurately model the DC, the nodes already computed are kept in the new sparse grid, thus reducing the number of evaluation of the DC coupling coefficient. (See Appendix 5.B for additional details on the Smolyak algorithm.) As a result, only 65 collocation points (Fig. 4) are required to build the desired SC model, and the values of the coupling coefficient at the interpolation nodes are computed using the FMM solver Fimmwave.
F. Directional Coupler Variability Analysis
Finally, the variability analysis for the coupling coefficient of the DC under study is performed using an SC model depending on the pair of independent Gaussian random variables and the results obtained are validated through comparison with an MC analysis based on the Fimmwave FMM solver on the DC cross section for the couple of correlated random variables . To compare the performance of the two methods, the same set of 10000 samples for the pair of correlated random variables (see Fig. 5). The corresponding values for the independent random variables are used to estimate the device variability features.
The proposed method shows an excellent accuracy compared with the classical MC analysis, as shown in Table 1, Figs. 6 and 7. In particular, the mean and the standard deviation of the coupling coefficient obtained employing the two methods are reported in Table 1: the relative error in the estimation of the mean and the standard deviation is only and , respectively. Apart from stochastic moments, more complicated functions of the stochastic process under study can be estimated: the probability density and cumulative distribution function (CDF) of obtained utilizing the two methods considered are in excellent agreement, as shown in Fig. 7.
From the field coupling coefficient, we can also easily derive performance parameters of a DC such as 3 dB coupling length from8, the proposed method shows excellent modeling accuracy and a good estimation of the PDF and the CDF of the 3 dB coupling length .
Furthermore, as presented in Table 2, the SC method has dramatically saved computational cost. Note that the SC method took a two-step procedure to perform the same variability analysis. Initially, SC required 65 simulations to compute the coupling coefficient at the collocation points. Next, we used the SC model over 10000 samples of the independent random variables in the MC method. Hence, the total computational time of the SC method is 8 min and 59 s, which represents a speed-up of a factor with respect to the MC analysis performed in Fimmwave for the couple of correlated random variables , which required 21 h 53 min 14 s.
This paper has presented the application of a novel technique for the efficient variability analysis of photonic devices, such as DCs. It is based on the use of SC methods. Thanks to the flexibility in the choice of the interpolation schemes and the efficiency of sparse grid sampling to choose the collocation nodes for multiple dimension, the proposed approach is flexible and can be applied to study a broad range of photonic devices. The accuracy and efficiency of the proposed technique have been verified using comparison with the classic MC analysis for a pertinent numerical example, achieving a simulation speed-up of .
KARHUNEN–LOÉVE EXPANSION AND CORRELATED GAUSSIAN RANDOM VARIABLES
Let us assume that the correlation matrix for the random variables under study is symmetric and positive-definite. Then, and can be diagonalized asA1), Eq. (9) becomes 9). Furthermore, it is possible to express the joint PDF Eq. (A2) with respect to a vector of independent Gaussian random variable , with zero mean and variance equal to , as 10) can be obtained by combining Eqs. (A4) and (A5).
Let us express a stochastic process depending on one random variable by means of the Lagrange interpolation scheme as 2). The nodes can be chosen from a node distribution, which guarantees a good quality of the approximation (i.e., the extrema of the Chebyshev polynomials). Extending Eq. (A6) to the case of multiple random variables can be performed via tensor product, as has been shown in Section 2, and Eq. (A6) becomes A6) with respect to the random variable , and is the number of random parameters considered. The total number of nodes required to compute Eq. (A7) is given by the product of the nodes used for each random parameter, as shown in Eq. (4). Clearly, the required number of nodes grows quickly with respect to the number of parameters considered. Indeed, if only two nodes are used for each random variable, the total number of points required for a full-tensor product interpolation is .
The Smolyak algorithm allows us to build multidimensional interpolation functions based on a minimal number of nodes by expressing the desired interpolation as a linear combination of tensor products. In particular, the property of the 1D interpolation is conserved for higher dimensions. Indeed, the sparse interpolant given by the Smolyak algorithm is
Let us denote as the set of points utilized in the 1D function interpolation. According to Eq. (A8), the stochastic process must be computed at the nodes of the sparse grid given by
1. M. Streshinsky, R. Ding, Y. Liu, A. Novack, C. Galland, A.-J. Lim, G. Q. Lo, T. Baehr-Jones, and M. Hochberg, “The road to affordable, large-scale silicon photonics,” Opt. Photon. News 24(9), 32–39 (2013). [CrossRef]
2. W. Bogaerts, M. Fiers, and P. Dumon, “Design challenges in silicon photonics,” IEEE J. Sel. Top. Quantum Electron. 20, 1–8 (2014). [CrossRef]
3. L. Chrostowski, X. Wang, J. Flueckiger, Y. Wu, Y. Wang, and S. T. Fard, “Impact of fabrication non-uniformity on chip-scale silicon photonic integrated circuits,” in Optical Fiber Communication Conference (Optical Society of America, 2014), paper Th2A.37.
4. S. K. Selvaraja, W. Bogaerts, P. Dumon, D. Van Thourhout, and R. Baets, “Subnanometer linewidth uniformity in silicon nanophotonic waveguide devices using CMOS fabrication technology,” IEEE J. Sel. Top. Quantum Electron. 16, 316–324 (2010). [CrossRef]
5. G. S. Fishman, Monte Carlo: Concepts, Algorithms, and Applications (Springer-Verlag, 1996).
6. D. Spina, F. Ferranti, G. Antonini, T. Dhaene, and L. Knockaert, “Efficient variability analysis of electromagnetic systems via polynomial chaos and model order reduction,” IEEE Trans. Compon. Packag. Manuf. Technol. 4, 1038–1051 (2014). [CrossRef]
7. P. R. Johnston, “Defibrillation thresholds: a generalised polynomial chaos study,” in Conference in Computing in Cardiology, Cambridge, MA (Sept. 7–10, 2014).
8. K.-K. K. Kim and R. D. Braatz, “Generalized polynomial chaos expansion approaches to approximate stochastic receding horizon control with applications to probabilistic collision checking and avoidance,” in IEEE International Conference on Control Applications, Dubrovnik, Croatia (Oct. 3–5, 2012).
9. D. Cassano, F. Morichetti, and A. Melloni, “Statistical analysis of photonic integrated circuits via polynomial-chaos expansion,” in Advanced Photonics (Optical Society of America, 2013), paper JT3A.8.
10. T. Weng, Z. Zhang, Z. Su, Y. Marzouk, A. Melloni, and L. Daniel, “Uncertainty quantification of silicon photonic devices with correlated and non-Gaussian random parameters,” Opt. Express 23, 4242–4254 (2015). [CrossRef]
11. M. S. Eldred, “Recent advance in non-intrusive polynomial-chaos and stochastic collocation methods for uncertainty analysis and design,” in Proceedings of 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, California (May , 2009).
12. N. Agarwal and N. R. Aluru, “Weighted Smolyak algorithm for solution of stochastic differential equations on non-uniform probability measures,” Int. J. Numer. Meth. Eng. 85, 1365–1389 (2011). [CrossRef]
13. N. Agarwal and N. R. Aluru, “Stochastic analysis of electrostatic mems subjected to parameter variations,” J. Microelectromech. Syst. 18, 1454–1468 (2009). [CrossRef]
14. D. Xiu, “Fast numerical methods for stochastic computations: a review,” Commun. Comput. Phys. 5, 242–272 (2009).
15. W. A. Weiser and S. E. Zarantonello, “A note on piecewise linear and multilinear table interpolation in many dimensions,” Math. Comput. 50, 189–196 (1988). [CrossRef]
16. B. Ganapathysubramanian and N. Zabaras, “Sparse grid collocation schemes for stochastic natural convection problems,” J. Comput. Phys. 225, 652–685 (2007). [CrossRef]
17. V. Barthelmann, E. Novak, and K. Ritter, “High dimensional polynomial interpolation on sparse grid,” Adv. Comput. Math. 12, 273–288 (2000). [CrossRef]
18. E. Novak and K. Ritter, “High dimensional integration of smooth functions over cubes,” Numer. Math. 75, 79–97 (1996). [CrossRef]
19. E. Novak and K. Ritter, “Simple cubature formulas with high polynomial exactness,” Constr. Approx. 15, 499–522 (1999). [CrossRef]
20. A. Der Kiureghian and P. L. Liu, “Structural reliability under incomplete probability information,” J. Eng. Mech. 112, 85–104 (1986). [CrossRef]
21. M. Loeve, Probability Theory, 4th ed. (Springer-Verlag, 1977).
22. W. Bogaerts, D. Taillaert, B. Luyssaert, P. Dumon, J. Van Campenhout, P. Bienstman, D. Van Thourhout, R. Baets, V. Wiaux, and S. Beckx, “Basic structures for photonic integrated circuits in silicon-on-insulator,” Opt. Express 12, 1583–1591 (2004). [CrossRef]
23. L. Chrostowski and M. Hochberg, Silicon Photonics Design: From Devices to Systems (Cambridge University, 2015).