## Abstract

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

## 1. INTRODUCTION

Integrated photonics, and in particular, silicon photonics, is rapidly enabling complex photonic functions on a chip [1]. 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 [4]. 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 [5] 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 [11].

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 $2\times 2$ 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 [12] (also defined as collocation points). Indeed, a stochastic process $\mathbf{Y}(\mathit{\xi})$ can be expressed as

where $\mathit{\xi}$ denotes the $N$ stochastic parameters and ${\{{L}_{i}(\mathit{\xi})\}}_{n=1}^{N}$ represents the interpolation basis functions.For a photonic device, the process $\mathbf{Y}$ could correspond to the functional parameters such as the waveguide propagating constants and the coupling coefficients in coupling devices. The stochastic variables ${\mathit{\xi}}_{i}$ 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 [15]). 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 ${L}_{i}$ in Eq. (1) for a 1D interpolation can be expressed as

where ${L}_{i}$ is equal to 1 for $\xi ={\xi}_{j}$ and is equal to 0 for $\xi ={\xi}_{i}$. Next, for interpolation in multiple dimensions, a tensor-product approach can be used, and Eq. (1) becomesThe 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 $\mathit{\xi}$ [11]. In the general case of correlated random variables, decorrelation can be obtained via a variable transformation, such as the Nataf transformation [20] or the Karhunen–Loéve expansion [21].

The stochastic moments (mean, variance, …) can be computed utilizing analytical formulas and then efficiently, once the analytical form of the interpolation functions ${\{{L}_{i}(\mathit{\xi})\}}_{n=1}^{N}$ has been decided. For example, if the random variables $\mathit{\xi}$ are defined in the sample space $\mathrm{\Omega}$, the mean of $\mathbf{Y}(\mathit{\xi})$ is defined as

*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 $\mathit{\xi}$ on the variations of the stochastic process considered $\mathit{Y}$ (dynamic stochastic processes require a higher number of collocation points);
- • the interpolation scheme ${L}_{i}$ 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 [16] 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. [22]. 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 (${n}_{\text{eff}}$) 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 $2\times 2$ 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 $z$. The power coupling $K(z)$ in a DC can be expressed as [23]

The power coupling consists of two parts: the contribution of the straight waveguide $\kappa z$ and of the bend section ${\kappa}_{0}$. The bend part will introduce an initial phase in the coupling term. Usually, the bend contribution is fairly small, and the power coupling is sinusoidally varying with the waveguide coupling section length $z$.The rate of coupling is defined as the field coupling coefficient $\kappa $, 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 $w$, thickness $t$, and gap $g$ (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 $g$ and $2\times $ the half-waveguide width $w$ is constant, as shown in Fig. 1. Therefore, in our example, we can describe the full geometry of the DC with only two parameters: $w$ and $t$.

In this study, we will use the SC technique to find out how geometry variability influences the DC performance, namely, the coupling coefficient $k$. Indeed, due to the fabrication variations, the fabricated linewidth $w$, thickness $t$, and gap $g$ 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 $w$ and thickness $t$ 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 $t$ and $w$ adhere to a Gaussian distribution.

#### B. Simulation Setup

According to the theory of supermodes, we can write the coupling coefficient $\kappa $ as [23]

where ${n}_{\mathrm{eff}\_o}$ and ${n}_{\mathrm{eff}\_e}$ are the effective index of asymmetrical and symmetrical supermodes in DC. For our silicon photonics devices, we assume the wavelength to $\lambda =1.55\text{\hspace{0.17em}}\mathrm{\mu m}$. Next, the nominal value of the width and thickness are ${w}_{0}=450\text{\hspace{0.17em}}\mathrm{nm}$ and ${t}_{0}=220\text{\hspace{0.17em}}\mathrm{nm}$, respectively, while we fix the sum of width $w$ and gap $g$ at 650 nm.To calculate $\kappa $ of a given geometry, we define the DC structure accordingly and simulate ${n}_{\mathrm{eff}\_o}$ and ${n}_{\mathrm{eff}\_e}$ 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 $K$ of a DC as a stochastic process depending on two correlated random variables with Gaussian PDFs: the width $w$ and thickness $t$. Hence, the joint PDF of the two random variables considered is defined as

To validate the robustness of the proposed method, ${\sigma}_{w}$ and ${\sigma}_{t}$ are chosen equal to 2% and the correlation coefficient $\rho =0.9$, 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 $\mathit{\eta}$, 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 [21], it is possible to express the vector of correlated Gaussian random variables $\mathit{\eta}$ in the vector of uncorrelated Gaussian random variables with zero mean and unit variance $\mathit{\xi}={[{\xi}_{1},{\xi}_{2}]}^{T}$ as

where $\mathbf{E}$ and $\mathbf{V}$ are the diagonal matrix of the eigenvalues and the full matrix of the eigenvectors of the covariance matrix $\mathbf{C}$, respectively. Because uncorrelated Gaussian random variables are also independent, we have now expressed the coupling coefficient as a stochastic process, which depends on the pair of independent Gaussian random variables $({\xi}_{1},{\xi}_{2})$. An accurate description of the Karhunen–Loéve expansion for Gaussian random variables is given in Appendix 5.A.#### 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: ${\xi}_{1}$ and ${\xi}_{2}$. In this example, the Clenshaw–Curtis rule is adopted [17]: 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 [17]. 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 $({\xi}_{1},{\xi}_{2})$ 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 $(w,t)$. To compare the performance of the two methods, the same set of 10000 samples for the pair of correlated random variables $(w,t)$ (see Fig. 5). The corresponding values for the independent random variables $({\xi}_{1},{\xi}_{2})$ 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 $9.0\times {10}^{-5}$ and $5.6\times {10}^{-3}$, 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 $\kappa $ 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 from

As shown in Fig. 8, the proposed method shows excellent modeling accuracy and a good estimation of the PDF and the CDF of the 3 dB coupling length ${l}_{3\text{\hspace{0.17em}}\mathrm{dB}}$.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 $146\times $ with respect to the MC analysis performed in Fimmwave for the couple of correlated random variables $(w,t)$, which required 21 h 53 min 14 s.

## 4. CONCLUSION

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 $146\times $.

## Appendix A

## KARHUNEN–LOÉVE EXPANSION AND CORRELATED GAUSSIAN RANDOM VARIABLES

Let us assume that the correlation matrix ${\mathbf{C}}^{N\times N}$ for the random variables $\mathit{\eta}$ under study is symmetric and positive-definite. Then, $\mathbf{C}$ and can be diagonalized as

Thanks to Eq. (A1), Eq. (9) becomes## SMOLYAK ALGORITHM

Let us express a stochastic process $\mathbf{Y}$ depending on one random variable $\xi $ by means of the Lagrange interpolation scheme as [16]

where ${L}_{i}$ is given by Eq. (2). The $Q$ 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) becomesThe 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 ${\mathit{A}}_{q,N}$ given by the Smolyak algorithm is

Let us denote $\mathrm{\Theta}$ as the set of points utilized in the 1D function interpolation. According to Eq. (A8), the stochastic process $\mathbf{Y}$ must be computed at the nodes of the sparse grid ${\mathbf{H}}_{q,N}$ given by

## REFERENCES

**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).