Stochastic process design kits for photonic circuits based on polynomial chaos augmented macro-modelling

Fabrication tolerances can significantly degrade the performance of fabricated photonic circuits and process yield. It is essential to include these stochastic uncertainties in the design phase in order to predict the statistical behaviour of a device before the final fabrication. This paper presents a method to build a novel class of stochastic-based building blocks for the preparation of Process Design Kits for the analysis and design of photonic circuits. The proposed design kits directly store the information on the stochastic behaviour of each building block in the form of a generalized-polynomial-chaos-based augmented macro-model obtained by properly exploiting stochastic collocation and Galerkin methods. Using these macro-models, only a single deterministic simulation is required to compute the stochastic moments of any arbitrary photonic circuit, without the need of running a large number of time-consuming circuit simulations thereby dramatically improving simulation e ciency. The e ectiveness of the proposed approach is verified by means of classical photonic circuit examples with multiple uncertain variables. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (250.5300) Photonic integrated circuits ; (000.5490) Probability theory, stochastic processes, and statistics (230.7370) Waveguides. References and links 1. X. Chen, M. Mohamed, Z. Li, L. Shang, and A. R. Mickelson, “Process variation in silicon photonic devices,” Appl. Opt. 52, 7638–7647 (2013). 2. D. Melati, F. Morichetti, A. Canciamilla, D. Roncelli, F. M. Soares, A. Bakker, and A. Melloni, “Validation of the building-block-based approach for the design of photonic integrated circuits,” J. Lightwave Technol. 30, 3610–3616 (2012). 3. Z. Lu, J. Jhoja, J. Klein, X. Wang, A. Liu, J. Flueckiger, J. Pond, and L. Chrostowski, “Performance prediction for silicon photonics integrated circuits with layout-dependent correlated manufacturing variability,” Opt. Express 25, 9712–9733 (2017). 4. A. Waqas, D. Melati, and A. Melloni, “Sensitivity analysis and uncertainty mitigation of photonic integrated circuits,” J. Lightwave Technol. 35, 3713–3721 (2017). 5. D. Cassano, F. Morichetti, and A. Melloni, “Statistical analysis of photonic integrated circuits via polynomial-chaos expansion,” in “Signal Processing in Photonic Communications,” (Optical Society of America, 2013), pp. JT3A–8. 6. A. Waqas, D. Melati, and A. Melloni, “Stochastic simulation and sensitivity analysis of photonic circuit through morris and sobol method,” in “Optical Fiber Communications Conference and Exhibition (OFC), 2017,” (IEEE, 2017), pp. 1–3. 7. M. Villegas, F. Augustin, A. Gilg, A. Hmaidi, and U. Wever, “Application of the polynomial chaos expansion to the simulation of chemical reactors with uncertainties,” Math. Computers Simulation 82, 805–817 (2012). 8. P. Sochala and O. Le Maître, “Polynomial chaos expansion for subsurface flows with uncertain soil parameters,” Adv. Water Resources 62, 139–154 (2013). 9. T.-W. 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). 10. Y. Xing, D. Spina, A. Li, T. Dhaene, and W. Bogaerts, “Stochastic collocation for device-level variability analysis in integrated photonics,” Photonics Res. 4, 93–100 (2016). Vol. 26, No. 5 | 5 Mar 2018 | OPTICS EXPRESS 5894 #313634 https://doi.org/10.1364/OE.26.005894 Journal © 2018 Received 21 Nov 2017; revised 18 Jan 2018; accepted 18 Jan 2018; published 26 Feb 2018 11. T.-W. Weng, D. Melati, A. Melloni, and L. Daniel, “Stochastic simulation and robust design optimization of integrated photonic filters,” Nanophotonics 6, 299–308 (2017). 12. R. Ghanem and P. D. Spanos, “A stochastic galerkin expansion for nonlinear random vibration analysis,” Probabilistic Engineering Mechanics 8, 255–264 (1993). 13. Z. Zhang, T. A. El-Moselhy, I. M. Elfadel, and L. Daniel, “Stochastic testing method for transistor-level uncertainty quantification based on generalized polynomial chaos,” IEEE Trans. Computer-Aided Design Integrated Circuits Systems 32, 1533–1545 (2013). 14. D. Xiu and G. E. Karniadakis, “Modeling uncertainty in flow simulations via generalized polynomial chaos,” J. Comput. Phys. 187, 137–167 (2003). 15. “Aspic from Filarate srl,” http://www.aspicdesign.com. Accessed: 30-10-2017. 16. M. S. Eldred, “Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design,” AIAA Paper 2274, 37 (2009). 17. D. Xiu and G. E. Karniadakis, “The Wiener-Askey polynomial chaos for stochastic di erential equations,” SIAM J. Sci. Comput. 24, 619–644 (2002). 18. D. Xiu, Numerical methods for stochastic computations: a spectral method approach (Princeton University Press, 2010). 19. D. Spina, F. Ferranti, T. Dhaene, L. Knockaert, G. Antonini, and D. V. Ginste, “Variability analysis of multiport systems via polynomial-chaos expansion,” IEEE Trans. Microwave Theory Techniques 60, 2329–2338 (2012). 20. J. B. Preibisch, P. Triverio, and C. Schuster, “Design space exploration for printed circuit board vias using polynomial chaos expansion,” in “Electromagnetic Compatibility (EMC), 2016 IEEE International Symposium on,” (IEEE, 2016), pp. 812–817. 21. J. B. Preibisch, P. Triverio, and C. Schuster, “E cient stochastic transmission line modeling using polynomial chaos expansion with multiple variables,” in “Numerical Electromagnetic and Multiphysics Modeling and Optimization (NEMO), 2015 IEEE MTT-S International Conference on,” (IEEE, 2015), pp. 1–4. 22. P. Manfredi, I. S. Stievano, and F. G. Canavero, “Parameters variability e ects on microstrip interconnects via hermite polynomial chaos,” in “Proc. of the 19th Conference on Electrical Performance of Electronic Packaging and Systems,” (2010), pp. 149–152. 23. P. Manfredi, D. V. Ginste, D. De Zutter, and F. G. Canavero, “Uncertainty assessment of lossy and dispersive lines in spice-type environments,” IEEE Trans. Components Packaging Manufacturing Technol. 3, 1252–1258 (2013). 24. D. Spina, T. Dhaene, L. Knockaert, and G. Antonini, “Polynomial chaos-based macromodeling of general linear multiport systems for time-domain analysis,” IEEE Trans. Microwave Theory Techniques 65, 1422–1433 (2017). 25. Y. Ye, D. Spina, P. Manfredi, D. V. Ginste, and T. Dhaene, “A comprehensive and modular stochastic modeling framework for the variability-aware assessment of signal integrity in high-speed links,” IEEE Trans. Electromagnetic Compatibility 60, 459–467 (2018). 26. A. Papoulis, “Random variables and stochastic processes,” (McGraw-Hill, 1985). 27. X. Leijtens, P. Le Lourec, and M. Smit, “S-matrix oriented cad-tool for simulating complex integrated optical circuits,” IEEE J. Sel. Top. Quantum Electron. 2, 257–262 (1996). 28. C. K. Madsen and J. H. Zhao, Optical Filter Design and Analysis: A Signal Processing Approach Optical Filter Design and Analysis: A Signal Processing Approach (Wiley Online Library, 1999). 29. A. Melloni and M. Martinelli, “Synthesis of direct-coupled-resonators bandpass filters for wdm systems,” J. Lightwave Technol. 20, 296–303 (2002).


Introduction
Recent advances in photonic integration are making possible the implementation of complex photonic circuits combining many functions on a single chip, significant production volumes and reduced fabrication costs [1]. Boosted by these progresses, Process Design Kits (PDKs) have recently emerged in photonics as a powerful tool for the analysis and design of complex circuits and as useful framework for in-depth exploitation of the potential of photonic for the large-scale integration. A PDK contains a large amount of foundry-specific information that can be used to design, simulate and fabricate circuits exploiting a given fabrication technology. Among these information, a PDK o ers to the users a library of the building blocks made available by the foundry, incorporating in particular a deterministic mathematical macro-model of the behaviour of each device, generally in the form of either a transmission or a scattering matrix [2]. In this scenario, a major issue is represented by fabrication tolerances, which can be particularly detrimental in large and complex circuits. Statistical variations in the fabrication processes of a device, such as waveguide width or height, proper gap opening or changes in material composition can have a strong e ect on its functionality and ultimately a ect the entire circuit performance [3][4][5][6]. The possibility to include information on the expected variability in each building block and the availability of e cient computational strategies to predict the statistical behaviour of the final circuit are hence of primary importance.
Monte Carlo is a robust, accurate, and easy-to-implement approach, and is often considered as the standard for stochastic analysis. However, its high computational cost can hamper its application to the analysis of even relatively simple circuits because a large number of simulation runs (generally in the order of 10 4 10 5 ) is required to obtain reliable results. For this reason, the generalised polynomial chaos approach has been introduced in several application fields as an e cient alternative to the classical Monte Carlo method [7,8], and recently it has been proposed also for the variability analysis of photonic devices [9][10][11]. The generalised polynomial chaos technique allows approximating the dependence of the simulation output on the stochastic input parameters with a set of orthonormal polynomials. The computation of the coe cients of the polynomial basis functions is the core task of this approach and can be done either through intrusive methods (e.g. stochastic Galerkin [12] and stochastic testing [13]), that require modifying the internal code of an existing deterministic solver, or non-intrusive sample-based methods (such as stochastic collocation [14]), that use the deterministic solvers as black boxes. The few Generalised Polynomial Chaos (gPC) implementations described in the literature for photonics applications are based on nonintrusive methods. Figure 1(a) schematically reports the common approach for stochastic analyses at circuit level. The building blocks available in the process design kit are exploited by the designer to build a circuit model [15] depending on a set of random parameters. During simulations, the value of these parameters is sampled according to their probability density function (e.g. Monte Carlo simulations) and the deterministic circuit simulator is used iteratively to compute the response of the designed circuit for each sample exploiting the deterministic model (e.g. scattering matrices) included in each building block. However, although the required number of simulations depends on the specific problem, it is generally smaller than few hundreds, as opposed to the case of direct Monte Carlo analysis. This small pool of initial Monte Carlo simulations is used to approximate the stochastic behaviour of the circuit through a generalized polynomial chaos expansion. This approach is obviously circuit-specific and has to be repeated each time the circuit parameters or its layout change.
In this work, we propose a Building-Block-based gPC method (BB-gPC) in which we exploit a generalized polynomial chaos expansion approach to realize a completely novel class of device models to be used within photonic PDKs. These stochastic models inherently convey stochastic information, which allows performing statistical analyses without any repeated simulation end enables an unprecedented simulation e ciency compared to classical gPC implementations. Stochastic collocation and stochastic Galerkin are combined to build augmented macro-models of each building block, directly embedding stochastic information in the form of gPC coe cients. The BB models, in the form of transmission or scattering matrices, are circuit independent and can be stored and replace the original deterministic macro-model of the building blocks in the process design kit [see Fig.1(b)]. The new matrices can hence be combined according to the building blocks connections to derive with a single run of the deterministic circuit simulator the stochastic behaviour of any circuit. The stochastic properties of the BB can reflect for example the foundry technological process and they are embedded in the PDK and should not be recalculated for every circuit.
The paper is structured as follows. A background overview of the generalised polynomial chaos approximation properties is given in Section 2, while the proposed method "BB-gPC" to build stochastic process design kits for photonic circuits is described in Section 3. The validation of the BB-gPC technique is performed in Section 4, assessing the accuracy and speed-up compared to direct Monte Carlo and classical gPC by means of two classical photonic circuit examples with multiple uncertain input variables. Conclusions and future perspective are summarized in Section 5.

Polynomial chaos background
The generalised polynomial chaos expansion is a well-known technique to describe finite-variance stochastic process Y depending on a vector of normalised random variables AE ⇠ as a summation of basis functions i ( AE ⇠) with suitable coe cients y i [16]: In this expression, AE ⇠ = [⇠ 1 , ⇠ 2 , · · · , ⇠ N ] are independent standard normal random variables, N is the number of random variables and i are orthogonal polynomials with respect to the probability measure W( AE ⇠), where ij is the Kronecker delta and a i are positive numbers. The construction of the gPC expansion in Eq. (1) entails the following three-step process. The first step is to determine the orthogonal polynomials. The second step is to truncate the series to a finite order for computational purposes, and the last step is to compute the gPC coe cients y i . Since the random variables AE ⇠ are assumed independent, the corresponding basis functions i ( AE ⇠) are computed as the product of the orthogonal polynomials corresponding to each individual random variable ⇠ i . It is important to note that for random variables with a specific distribution (i.e., Gaussian, Uniform, Beta, etc.), the polynomials tabled in the Wiener-Askey scheme can be used as basis functions, yielding an optimal (exponential) convergence rate [17]. For example, for the variables with Gaussian probability density function (PDF), the basis functions are the Hermite polynomials, while for the uniform PDF the basis functions are the Legendre polynomials. Considering all N-dimensional polynomials ( AE ⇠) of order P (total-order expansion), the expansion (1) can be truncated to M + 1 terms with the number of gPC basis function given by [18] M After the determination of the basis functions and truncation of the series in Eq. (3), the M + 1 scalar coe cients y i can be computed, for example, with a nonintrusive approach based on linear regression [16]. An initial small pool of Monte Carlo simulations is generated for the process Y and the coe cients y i are calculated in order for Eq. (3) to match the known solutions in a least-mean-square sense. In this work, the fabrication uncertainty is represented as independent and Gaussian-distributed random variables. In case linear correlation exist between variables, the problem could be solved by decorrelation techniques and will be treated in future works. Also, we consider as the stochastic process under study the complex frequency response (transfer function) of photonic devices and circuits, which can be handled by considering Eq. (3) for each wavelength and calculating the coe cients y i as wavelength dependent. Alternatively, another method to handle the frequency dependence is a rational approximation of the frequency response through the vector fitting algorithm to obtain a state-space model of the transfer function that can then be approximated through Eq. (3) [5,19].

Augmented macro-model via polynomial chaos expansion
The polynomial chaos formalism described in the previous section is exploited here to create an augmented-macro model of a building block embedding its stochastic behaviour. The formalism presented in this section is based on the concept of augmented macro-modelling which has recently been introduced in electronics and microwave for the stochastic analysis of printed circuit board [20], transmission line [19,21], micro-strip interconnects [22][23][24] and high-speed links [25]. The description provided is based on scattering matrix formalism and is only valid for linear and passive photonic devices. However, the same approach can be immediately extended also to transfer matrices. For a given building block, the scattering matrix represents the relation between the complex amplitude of the incoming fields waves and the complex amplitude of the outgoing waves. In this formalism, the behaviour of a building block with n p ports is expressed as where vectors AE a 2 C n p ⇥1 and AE b 2 C n p ⇥1 are the complex input and output waves at all the n p ports and AE S 2 C n p ⇥n p is the scattering matrix containing all the scattering parameters of the building block. In the case of stochastic behaviour, the scattering matrices of the BB and the complex wave amplitudes depends on the set of random variables AE ⇠. Using the gPC expansion (3), the wave relation (5) becomes where AE a i , AE b i , AE S i are now the polynomial chaos coe cients of the input and output wave amplitudes and of the building block scattering parameters for each port, respectively. A non-intrusive approach can be used to obtain the gPC coe cient AE S i of a given building block. Similarly to the method described in section 2, in this work we exploit Monte Carlo simulations sampling an initial set of Q discrete samples of the normalised random variables AE The gPC coe cients are calculated solving the least square system where the j th row of the matrix contains the multivariate polynomial basis evaluated at AE ⇠ j and the j th column of the matrix R represents the corresponding complex response of the building block. In order to build the augmented scattering matrix of the building block, the Galerkin method [16,19] is used. The projection of Eq. (6) on the p th gPC basis function p gives where the elements AE b k p represent the p th gPC coe cients at the k th port (k = 1 . . . n p ) for the output wave. The factors are real scalar numbers, which can be computed upfront solving multi-dimensional integrals analytically only once and then stored [26]. It is worth mentioning that most of the integrals are zero and therefore the resulting matrix is highly sparse, limiting the computational cost. Repeating this operation for all the M + 1 gPC basis functions leads to where the vectors AE a PC 2 C (M+1)n p ⇥1 and AE b PC 2 C (M+1)n p ⇥1 contain the gPC coe cients of the input and output waves respectively, while the augmented matrix AE S PC 2 C (M+1)n p ⇥(M+1)n p of the building block is a weighted combination of gPC coe cients. As an example, Fig. 2(a) shows a standard photonic circuit in which the port A of BB 1 (ring) is connected with port B of BB 2 (directional coupler). Figure 2(b) shows the augmented representation of the standard circuit in which ports A 0 , A 1 and A 2 are connected with B 0 , B 1 and B 2 , respectively. Note that the augmented circuit has no physical meaning. The augmented BBs are larger as each port is scaled by the factor (M+1), however the corresponding matrices are highly sparse with a large number of null entries. The additional overhead in handling the sparse augmented matrices, at least for few tens of parameters, is negligible compared to the time required to run a Monte Carlo analysis. In order to clarify the procedure adopted to compute AE S PC , let us assume that the building block under study depends only on one random variable and the gPC expansion of its scattering parameters is formed by two basis functions. In this case Eq. (6) becomes For the sake of clarity, the dependence on the random vector AE ⇠ is omitted. Due to the orthonormality relation in Eq. (2), the projection of Eq. (10) onto the first basis function 0 yields Further, projecting Eq. (10) onto the second basis function 1 gives Upon calculation of the scalar products in Eqs. (11) and (12), the augmented matrix relating the expansion coe cients of the input to the coe cients of the output is derived as with Hence, Eq. (13) describes the "augmented model" of the building block which depends on the considered random variables. It is important to mention that AE S PC is still a unitary and loss-less scattering matrix and the use of orthonormal gPC basis functions preserves the symmetry of AE S PC , as can be seen from Eq. (14). Next, the procedure can be repeated for all the needed frequency points in a given frequency range to obtain frequency-dependent augmented models AE S PC ( l ) for l = 1, · · · , L. Note that the expansion is calculated independently at each wavelength but starting from the initial set of Q Monte-Carlo samples that are in the form of wavelength-dependent transfer function thus including the dependence on wavelength in the augmented model. Alternatively, a frequency model based for example on vector fitting can be exploited. Moreover, if the original matrix AE S depends also on a set of non-random design parameters, the gPC coe cients can be calculated for di erent values of these parameters. The coe cients of the gPC approximation are then fitted in the allowed range of the design parameters, to obtain a parametric augmented model.
The library of BBs modeled by frequency-dependent augmented matrices in form of Eq. (6) realize a Stochastic Process Design Kit. Using standard procedures [2,27], the BBs augmented scattering matrices can be combined according to the circuit layout with a single simulation run of available deterministic circuit simulators and the stochastic behaviour of any arbitrary circuit is obtained.

Numerical Results
In this section, the proposed approach is applied to compute the augmented models of two di erent building blocks that are hence exploited to analyze the stochastic behaviour of two di erent photonic circuits, an unbalanced second-order Mach-Zehnder filter and a fifth-order coupled ring resonator filter. In the following we consider standard silicon-on-insulator technology. In order to demonstrate the potential of the proposed approach, circuits are assembled combining the same standard elementary building blocks: the directional coupler and the single waveguide. The augmented macro-models of these two building blocks are calculated once and stored as described in Section 3. The complex frequency response (transfer function) of the photonic filters is considered as the stochastic process under investigation. The fabrication uncertainty is represented as independent and Gaussian-distributed random variables and the corresponding basis functions are products of Hermite polynomials.
In order to calculate the augmented macro-models, we computed the following steps as described in Section 3: The obtained wavelength-dependent gPC coe cients allow to calculate the stochastic moments, such as mean and variance, simply as [14] µ( ) = AE S 00 PC ( ), Stochastic functions such as the probability density function (PDF) and the cumulative density function (CDF) can be computed following standard analytical formulae or by applying Monte Carlo sampling on the gPC approximation of the circuit under study, described by the scattering coe cients In order to validate the reliability of the proposed technique, two classical and well known filters, a second-order cascaded Mach-Zehnder filter and a fifth-order coupled ring resonator filter are considered. We compare the obtained results by BB-gPC with those obtained by classical gPC method as described in Section 2 and by Monte Carlo analysis. The total number of samples chosen for the Monte Carlo analysis represents a trade-o between accuracy and associated computational cost in estimating not only simple stochastic moments, such as mean and standard deviation, but also the probability density function. The simulations are performed with the commercial circuit simulator ASPIC [15] on a laptop with an Intel Core i5 processor at 2.4 GHz and 8-GB RAM.

Example A: Unbalanced Mach-Zehnder filter
The second-order Mach-Zehnder filter is assembled using two standard elementary building blocks, BB W and BB K , as schematically shown in Fig. 3(a). The nominal design of the filter was obtained with the synthesis technique described in reference [28]. The Mach-Zehnder filter has a nominal 3-dB bandwidth of BW 0 = 62 GHz and the corresponding coupling coe cients for the three directional couplers are K 1 =K 3 =0.865 and K 2 =0.532, with in-band isolation larger than 20 dB. For the directional coupler building block, the considered nominal gap distance is g i =0.3 µm, leading to coupling lengths of L c1 =L c3 =23.45 µm, and L c2 =16.09 µm. For the waveguide building block the nominal waveguide width is W i =408 nm and thickness is 220 nm, corresponding to an e ective index and group index of about 2.23 and 4.402, respectively. Both have the same unbalance lengths of 680.8 µm, corresponding to a free spectral range of 100 GHz. Figure 3(b) shows the ideal transfer function of the Mach-Zehnder filter at the bar (bold black) and the cross (bold red) ports respectively. These spectral responses do not take into account the unavoidable process variations a ecting real fabricated circuits. In the stochastic analysis, the waveguide width uncertainty ( W i ) and the gap uncertainty ( g i ) are assumed as independent normal distributed random variables with standard deviation W =10 pm and g =5 nm, respectively. Figure 3 Fig. 1(a)) at each wavelength.

Figures 4(a) and 4(b)
show the mean and the standard deviation of the transmission at the bar (black) and the cross (red) ports of the filter. Blue dash-dot lines represent the result obtained by BB-gPC while full line and circles shows the result obtained using Monte Carlo and gPC analysis, respectively. The proposed BB-gPC technique is in excellent agreement compared with the classical MC and gPC analysis in computing mean and standard deviation of the circuit. Figure 5(a) and 5(b) show the probability density function of the intensity transfer function at the bar (black) and the cross port (red) of the filter at the centre wavelength 1.55025 µm, respectively. Blue dash-dot line represents BB-gPC while full line and circles shows the results of classical MC and gPC respectively. Note that for both ports, the probability density function is asymmetric. The cross port of the filter su ers from a broader deviation of the intensity, while the bar port is more robust for the same process uncertainties. The obtained gPC approximation (17) can be used to perform the analysis of other desired quantities of interest, for example the 3-dB bandwidth of the filter. The transfer functions obtained with MC analysis on the original circuit and on both classical gPC and BB-gPC approximations are used to calculate the PDF of the 3-dB bandwidth of the filter at the bar port as shown in Fig. 5(c). The nominal bandwidth BW 0 is marked with a dashed line for reference. The bandwidth PDF is asymmetric. The obtained results are also exploited to estimate the yield that can be achieved by the process (percentage of fabricated devices respecting a certain performance constraint). For instance, considering the bandwidth constraint (BW 0 ± 2 GHz), the achievable device yield for the considered uncertainties is equal to 89 %. The results obtained with all the three methods are in good agreement demonstrating that the BB-gPC provides reliable results not only for simple stochastic moments such as mean and variance but also for complex stochastic moments as the probability density function of any quantity of interest. The computational time of this test case is summarized in Table 1. Advantageously, the proposed method required only 5.6 seconds to build augmented model, to retrieve gPC approximation of the entire circuit and to run 10 4 MC samples analysis. In addition, it requires 3 minutes and 9 seconds to compute gPC model using initial 30 samples for each BB in the form of Eq. (6). However, this computation is required to be performed only once to create the library (PDK), which is then used to obtain the statistical information of any circuit layout.

Example B: Coupled ring-resonators filter
In this example we exploit the same stochastic building blocks models BB K and BB W described and calculated in Section 4.1 to analyze a fifth-order coupled ring filter schematically shown in Fig.6(a). The filter is designed according to the synthesis technique described in reference [29]. With a nominal bandwidth BW 0 =28.5 GHz and FSR=100 GHz, the resulting coupling coe cients for the six directional couplers are K 1 =K 6 =0.82, K 2 =K 5 =0.37 and K 3 =K 4 =0. 22 providing an in-band isolation larger than 25 dB. Also in this case, silicon photonics technology is considered, with channel waveguide cross section equal to 408⇥220 nm, corresponding to an e ective index and group index of about 2.23 and 4.402, respectively, as described in Section 4.  as Section 4.1, leading to coupling lengths of L c1 =L c6 =22.29 µm, L c2 =L c5 =12.90 µm and L c3 =L c4 =9.6323 µm. Figure 6(b) shows the ideal transfer function of the ring filter at the drop (bold black) and the through (bold red) ports respectively. In this example, BB K is considered deterministic and only BB W is a ected by stochastic variations, therefore the number of random variables is N = 5 and the augmented matrix size is 42 ⇥ 42 (M = 21 and n p = 2). Although in practice also the couplers should be considered stochastic, this example demonstrate that is it possible to combine stochastic BBs with deterministic BBs. The computed augmented model of BB K is still valid, but only the first coe cient of the gPC approximation in the augmented model is actually used, whereas the remaining coe cient are set equal to zero as the considered standard deviation for coupling gap is g =0. For BB W , the waveguide widths ( W i ) are assumed as independent standard Gaussian with standard deviation W =10 pm. The thin grey lines in Fig. 6(b) show the e ect of process uncertainties on the coupled ring resonator filter obtained by Monte Carlo analysis. Similar to Section 4.1, the BB K and BB W are combined to obtain the augmented matrix description of the whole circuit and the gPC approximation in the form of Eq. (17). Using BB-gPC, the circuit statistical behavior is computed and comparison with Monte Carlo and classical gPC analysis is performed. Results for the stochastic moments of the ring filter are shown in Figs. 7-9. In particular, Figs. 7(a) and 7(b) show the mean and standard deviation of the coupled ring resonator filter at the drop (black colour) and the through (red colour) ports respectively. Blue dash-dot lines represent the result obtained by BB-gPC, while full line and circles shows the result obtained using Monte Carlo and gPC analysis, respectively. Figs. 8(a) and 8(b) show the PDF of the intensity at the drop and through ports of the filter at the centre wavelength and Fig. 8(c) shows the PDF of 3-dB bandwidth obtained from direct Monte Carlo analysis (full line), classical gPC analysis (circles) and BB-gPC method (blue dash-dot line). The probability density functions of both drop and through ports of the ring filter are asymmetric. Also in this case, the intensity at the through port su ers form a broader deviation when compared to drop port of the filter for the same process uncertainties. In this example, to further demonstrate the potential of the BB-gPC, the statistical behaviour of the group delay at the drop port is  calculated and shown in Figs. 9(a) and 9(b). The inset in Fig. 9(a) shows the PDF at the centre wavelength. Figures 7-9, clearly show that the results of our proposed method are in excellent agreement with Monte Carlo and classical gPC technique. The results also shows that once the stochastic augmented macro-model of BBs are build and stored in the PDK, they can be used to perform the variability analysis of arbitrary photonic circuit by means of only one single simulation run.

Conclusion
In this paper, we presented a novel class of stochastic process design kits for the e cient variability analysis of photonic circuits. In the proposed framework generalised polynomial chaos expansion is exploited to build augmented macro-models of each building block directly embedding the description of its stochastic behaviour. The augmented macro-models are computed through the stochastic collocation and Galerkin methods and are part of the library (PDK). Using these augmented models, the stochastic moments of an arbitrary photonic circuit can be evaluated by a single deterministic simulation run, without the need for repeated simulations. The e ciency and flexibility of the framework have been illustrated using relevant numerical examples. Even if in this work we exploited synthetic simulation data to realize the stochastic process design kit, real experimental data can also be used. The proposed method can be helpful to analyse the stochastic variability and perform robust design optimization to obtain the circuit design which is more robust to fabrication process variation and achieve higher yield respecting certain performance constraints. Additionally, the BB-gPC method can be also be advantageous to perform sensitivity analysis, identifying the most critical parameters a ecting circuit performances.