Stochastic simulation and robust design optimization of integrated photonic filters

Tsui-Wei Weng; Daniele Melati; Andrea Melloni; Luca Daniel

doi:10.1515/nanoph-2016-0110

Artikel Open Access

Stochastic simulation and robust design optimization of integrated photonic filters

Tsui-Wei Weng , Daniele Melati , Andrea Melloni und Luca Daniel

Veröffentlicht/Copyright: 30. Juli 2016

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Manuskript einreichen Informationen für Autor*innen

Aus der Zeitschrift Nanophotonics Band 6 Heft 1

Abstract

Manufacturing variations are becoming an unavoidable issue in modern fabrication processes; therefore, it is crucial to be able to include stochastic uncertainties in the design phase. In this paper, integrated photonic coupled ring resonator filters are considered as an example of significant interest. The sparsity structure in photonic circuits is exploited to construct a sparse combined generalized polynomial chaos model, which is then used to analyze related statistics and perform robust design optimization. Simulation results show that the optimized circuits are more robust to fabrication process variations and achieve a reduction of 11%–35% in the mean square errors of the 3 dB bandwidth compared to unoptimized nominal designs.

Keywords: design optimization; generalized polynomial chaos; process variations; photonic integrated circuits; yield; coupled ring resonator filters

1 Introduction

Photonics is rapidly emerging as a mature and promising technology, and it is evolving from a pure research topic to a market-ready player, aiming at achieving large production volumes and small fabrication costs. Pushed by these motivations, process design kits (PDKs), circuit simulators, generic foundry approaches, and multiproject wafer runs are quickly changing the way that photonic circuits are conceived and designed [1], [2], [3], [4].

On the contrary, stochastic uncertainties related to fabrication variations, such as waveguide geometry deviation, gap opening issues, material composition fluctuations, and surface roughness, are unavoidable in production processes [5], [6], [7], [8], [9]. It is well known that such uncertainties can have a dramatic impact on the functionality of fabricated circuits [4], [10], [11], [12], [13], [14], [15]. To obtain a high-quality design of a photonic circuit (e.g. high yield or smaller performance variability), it is important to include such uncertainties during the early design stages. Hence, uncertainty quantification techniques become fundamental instruments to efficiently obtain the statistical information of the circuits as well as to achieve a high-quality design.

Monte Carlo is an approach commonly exploited to evaluate the impact of fabrication uncertainties on the functionality of the designed circuits [12]. Although effective, it suffers from a slow convergence rate and requires long computation time. Meanwhile, stochastic spectral methods have recently been regarded as a promising alternative for statistical analysis due to their fast convergence. The key idea is to approximate the output quantity of interest (e.g. the power consumed by the circuit or the bandwidth of a filter) with a set of orthonormal polynomial basis functions, known as generalized polynomial chaos (gPC) expansion. There are two classes of method to compute the coefficients of the basis functions, and each class has its own pros and cons. For intrusive methods (i.e. nonsampling methods) such as stochastic Galerkin [16] and stochastic testing [17], the computation cost is sometimes lower but it requires modifying the internal code of an existing deterministic solver. Conversely, nonintrusive methods (i.e. sample-based methods), including stochastic collocation [18] and least-squares regression techniques, use the deterministic solvers as a black box, which is often more convenient in practice. In addition, if the problem at hand happens to be inherently sparse, a sparse gPC model can be constructed by minimizing the ℓ₁ norm of the gPC coefficients [19], [20], [21] or using a tensor recovery model with low-rank and sparse constraints [22].

When performing design optimization based on gPC models, one approach is to use evolution algorithms such as genetic algorithm, which optimizes circuit performance under both uncertainties and design constraints with the obtained gPC surrogate models [23], [24], [25], [26]. The cost functions are usually expected values and/or variance (and/or their combinations) of the quantity of interests. However, that approach must regenerate gPC models for every design point inside the optimization loop because an analytical cost function is unavailable. On the contrary, a combined gPC model is proposed in [27], [28] to expand the cost function in terms of both process variation variables and design variables. In that approach, the cost function can be expressed in terms of design variables; therefore, gradient-based optimization algorithms can be applied. In particular, the cost function is a multivariate polynomial because the gPC bases are polynomials; hence, a global polynomial optimization solver, such as gloptipoly3 [29], can be employed. It is shown in [29] that the optimizer will provide certificates once a global optimum is found.

In this work, we exploit the sparsity structure present in photonic circuits and construct a sparse combined gPC model to analyze their related statistics and perform design optimization. Section 2 briefly reviews some background material on gPC and techniques for building a combined gPC model. The idea of constructing a sparse gPC model for design optimization is illustrated in Section 3, and a real-world photonic coupled ring resonator filter circuit example is simulated and demonstrated to prove the effectiveness of the method in Section 4. Finally, conclusions are summarized in Section 5.

2 Background review

In this section, a brief background review of gPC is first given, and the idea of using gPC models in design optimization is described afterwards.

2.1 gPC model

Let u(x→, ξ→) be a quantity of interest of a device/circuit/system under process variations, for example, the bandwidth of filter, the quality factor or resonant wavelength of a resonator, and the effective phase index of a waveguide. The vector x→ denotes a set of constant design variables, such as the gap of a directional coupler, the nominal radius of a ring resonator, or the nominal width of a waveguide, and ξ→ is a random vector describing fabrication process variations, such as the variations of the gap, the radius, and the width. Given the joint probability density function (pdf) of ξ→, the idea of gPC is to approximate u(x→, ξ→) by a set of specialized basis functions:

(1)u(x→, ξ→)≈∑α→Cα→(x→)Ψα→(ξ→)

where Ψα→(ξ→) is a multivariate polynomial and Cα→(x→) is the corresponding coefficient with nonnegative integer multi-index α→=(α1, ⋯, αd). Assuming all components of ξ→ are independent to each other, the joint pdf pξ→ is a product of one-dimensional pdf pξi, and Ψα→(ξ→) is a product of univariate polynomials ψαi(ξi). For some well-known distribution, such as uniform, Gaussian, β, and γ distribution, the corresponding bases are also known, and they are Legendre, Hermite, Jacobi, and Laguerre polynomial, respectively. For more details, readers are referred to [18].

If a total-order truncation of α→ with order p is used, α→ satisfies

α1+⋯+αd≤p

and there are total N=(d+p)!p!d! terms in (1). The multi-indices can be mapped into a set of integers from 1 to N, so (1) can be rewritten as

u(x→, ξ→)≈∑n=1Ncn(x→)Ψn(ξ→).

The coefficients Cα→ can be obtained either intrusively by stochastic Galerkin method or sample-based method such as stochastic collocation method and regression method.

2.2 Design optimization with gPC models

It is of high interest to design a device or circuit that can still perform well under process variations. In other words, the robustness of its performance is of primary concern, because fabrication variations are often unavoidable in reality. The goal is to optimize the quantity of interest (or a function of the quantity of interest) under uncertainties and design constraints. Because the quantity of interest is a random variable, it is reasonable to use its expectation or its associated function as the objective in the optimization problem. A common form of the design optimization problem is the following:

(2)minimizex→Eξ[u(x→, ξ→)] and/or Vξ[u(x→, ξ→)]subject toai≤xi≤bi, i=1, …, m,

where u(x→, ξ→) is the quantity of interest. The expectation and variance of u over ξ→ are known to be

(3)Eξ[u(x→, ξ→)]=C0→(x→)

and

(4)Vξ[u(x→, ξ→)]=∑α→≠0→Cα→2(x→),

assuming x→ is a constant vector and that the basis functions Ψα→(ξ→) are orthonormal.

If (3) and (4) are used to solve the optimization problem (2), then u(x→, ξ→) has to be reexpanded each time for each design point x→ inside the optimization loop, which is time consuming and only non-gradient-based optimization algorithm such as evolutionary genetic algorithm can be used to search for the optimum. However, genetic algorithm is a heuristic-based method, and it does not guarantee local and global optima, although it is often used as a tool to try to find a global optimum in practice. Alternatively, if u(x→, ξ→) can be expanded both on x→ and ξ→, then (3) and (4) will be analytical functions of x→; hence, gradient-based optimizers can be employed. Notice that the constraints of x_i in (2) are boxed constraints, which is equal to express x_i as uniformly distributed variables in the interval [a_i, b_i]. Therefore, rewrite u(x→, ξ→) as

u(x→, ξ→)≈∑α→Cα→Ψα→(x→, ξ→)=∑β→,γ→Cβ→,γ→Lβ→(x→)Φγ→(ξ→)

where Lβ→ is a multivariate Legendre polynomial with multi-index β→, Φ is a multivariate polynomial of ξ→ with multi-index γ→, and α→=(β→, γ→).u(x→, ξ→) is now referred to the combined gPC model. Using the combined gPC model, we have

Eξ[u(x→, ξ→)]=∫u(x→, ξ→)ρξ→dξ→=∫∑β→,γ→Cβ→,γ→Lβ→(x→)Φγ→(ξ→)ρξ→dξ→=∑β→,γ→Cβ→,γ→Lβ→(x→)∫Φγ→(ξ→)ρξ→dξ→

Because Φγ→(ξ→) are orthogonal and Φ0→(ξ→)=1, we have

∫Φγ→(ξ→)ρξ→dξ→={0, if γ→≠0→1, if γ→=0→

Therefore, we have

(5)Eξ[u(x→, ξ→)]=∑β→Cβ→,0→Lβ→(x→).

Similarly, using the orthogonality of basis functions

∫Φγ→(ξ→)Φγ→′(ξ→)ρξ→dξ→={0, if γ→≠γ→′1, if γ→=γ→′,

we have

(6)Vξ[u(x→, ξ→)]=Eξ[u2(x→, ξ→)]−(Eξ[u(x→, ξ→)])2=∑β→,β→′,γ→≠0→Cβ→,γ→Cβ→′,γ→Lβ→(x→)Lβ→′(x→).

3 Our proposed method

We propose to exploit both the idea of sparse gPC model [21] and the combined gPC model [27], [28] to construct a sparse combined gPC model for the quantity of interests in photonics applications. Specifically, our example circuit is a five-ring coupled resonator filter, and the quantity of interest is the 3 dB bandwidth. The combined gPC model is an efficient way to perform design optimization; however, if the number of parameters (i.e. x→ and ξ→) is more than 10 (there are a total of 17 parameters in our photonics example), using stochastic collocation to compute gPC coefficients is very costly. On the contrary, the coefficients of the sparse gPC model can be computed efficiently by solving an ℓ₁ minimization with the off-the-shelf ℓ₁ minimization solver, such as spgl1 [30]. With several hundreds of parameters, spgl1 can give a solution in minutes. In addition, the number of terms of the cost function in (2) is smaller by exploiting the sparsity of the gPC expansion, which can help saving some evaluation time during optimization process.

3.1 Sparse combined gPC model

Our sparse combined gPC model is constructed as follows:

BW(x→, ξ→)≈∑n=1NcnΨn(x→, ξ→),

where we compute the coefficients c_n by solving the following problem:

(7)min∥c→∥1, subject to ∥Ψc→−BW∗∥2≤δ,

where c→=(c1, ⋯, cN) is a vector of all coefficients, Ψ is the measurement matrix whose (i, j) entry is the jth basis with ith sample Ψj(x→(i), ξ→(i)), BW*=(BW⁽¹⁾, …, BW^(M)) is a vector of M output samples obtained from Monte Carlo simulation, and δ is a bound on the data error, which is usually estimated from cross-validation data.

3.2 Robust design optimization under process variations

In our example, we would like the bandwidth to be as robust as possible to the fabrication variations. In other words, we would like to minimize the expected mean square error (MSE) of the bandwidth with respect to the original designed bandwidth (which is also called nominal bandwidth BW₀ in this paper). The analytical expression of the expected MSE is

(8)MSE(x→)=Eξ[(BW(x→, ξ→)−BW0)2]=Vξ[BW(x→, ξ→)]+(Eξ[BW(x→, ξ→)]−BW0)2,

where the mean and variance are shown to be a multivariate polynomial in (5) and (6), respectively.

Thus, the robust design optimization that we will be solving is the following:

(9)minimizex→MSE(x→)subject toai≤xi≤bi, i=1, …, m,

Because the objective function MSE is a multivariate polynomial, we can obtain the global optimum of a polynomial optimization problem by solving generalized problems of moments [29], [31], which is an additional benefit of using the combined gPC model in the design optimization problem. Figure 1 summarizes the design flow of the proposed technique. Note that we use BW to denote the quantity of interest (3 dB bandwidth) here to be consistent with our example in the next section, but it can be any other quantity based on different applications.

Figure 1:

Proposed robust design optimization flow.

4 Example of integrated photonic filters

4.1 Benchmark description

To demonstrate the application of our approach, a fifth-order directly coupled ring resonator filter is used as a test case, and its schematic is shown in Figure 2A, where n_eff,i denotes the effective phase index of each ring and gi is the gap width of each directional coupler. Two “nominal designs” of the filter are obtained with standard synthesis techniques described in [32], [33].

Figure 2:

(A) The example circuit is a five-ring coupled resonator filter. (B) Transfer functions of the drop port (red bold solid line) and through port (blue bold dashed line) of the Chebyshev filter nominal design. Gray thin lines plot 100 Monte Carlo simulations of the transfer functions when fabrication variations exist on the effective index of each ring and gap width of each directional coupler. The variations of the effective phase indices and the gaps are assumed to be zero-mean Gaussian distribution with standard deviation σ_neff=10^-5 and σ_g=5×10^-3 μm, respectively.

The first design is a Chebyshev type I filter (equi-ripple bandpass response) with in-band isolation of 26 dB at the through port. The coupling coefficients of the directional couplers are

K1,0C=K6,0C=0.337,K2,0C=K5,0C=0.024,K3,0C=K4,0C=0.012.

The filter has a 3 dB passband bandwidth BW₀=25.6 GHz. The design and transfer functions of the filters were calculated with a commercial circuit simulator [34]. The transmission of Chebyshev nominal design is shown in Figure 2B, where the red bold solid line is for drop port and the blue dash line is for through port. The second design is a Butterworth filter (maximally flat passband response) with the same 3 dB bandwidth BW₀=25.6 GHz. The corresponding coupling coefficients are

K1,0B=K6,0B=0.490,K2,0B=K5,0B=0.042,K3,0B=K4,0B=0.013,

with an in-band isolation larger than 50 dB. For both designs, all the rings have the same length of 336.2 μm, the gaps g_i are 0.3 μm, and the effective phase indices n_eff,i are 2.23. The free spectral range is 400 GHz. Without loss of generality, dispersion and waveguide attenuation are neglected. Although these results already represent optimum filter designs, they do not take into account the unavoidable process variations affecting real fabricated circuits. Therefore, our goal is to include such variations during the design phase. To this purpose, the effective indices and gaps are written as the sum of nominal values (n_eff,0, g_0,i) and variations (Δn_eff,i, Δgi) respectively:

neff,i=neff,0+Δneff,i,gi=g0,i+Δgi.

The variables Δn_eff,i and Δg_i are denoted as process variation variables in the combined gPC model, and they are assumed to be Gaussian distributed with zero mean and standard deviation σ_neff and σ_g, respectively. On the contrary, g_0,i are denoted as design variables, which lie uniformly in an interval with lower bound g_lb and upper bound g_ub. All the process variation variables and design variables are assumed to be independent.

Note that, when we refer to “nominal design” in this paper, we mean that it is the design obtained without fabrication variations. Hence, the two nominal designs described earlier have n_eff,0=2.23, g_0,i=0.3 μm, Δn_eff,i=0, Δg_i=0. When a nominal design is exposed to process variations, we will explicitly describe it as “nominal design with fabrication variations”. Such variations could heavily affect the function of the circuit, as shown in Figure 2B. The thin gray lines in Figure 2B plot the effect of fabrication uncertainties on the Chebyshev nominal design considering Δn_eff,i, Δg_i≠0, σ_neff=10^-5, and σ_g=5×10^-3. It is clearly seen that the transfer function of the original nominal design (bold lines) is heavily distorted. The 3 dB passband as well as the in-band isolation have large fluctuation even under relatively small process variations. Hence, our goal is to find the best nominal design whose performance is the most robust to the process variations by exploiting the approach described in Section 2. Specifically, our target is to find the best gap widths g_0,i, which minimize the fluctuation of 3 dB bandwidth in the two nominal designs with process variations. The gaps values are constrained with some specified lower bound g_lb and upper bound g_ub.

4.2 Simulation results and discussion

In this section, we describe the simulation procedures and results in detail. The whole design flow is summarized in Figure 1, and all simulations are performed on an Intel i5-5200 CPU laptop with 8 GB of RAM.

The quantity of interest is the filter’s 3 dB bandwidth BW (x→, ξ→). The design variables are x→=g→0=(g0,1, ⋯, g0,6), and the (zero-mean Gaussian distributed) process variation variables are ξ→=(Δg1, ⋯, Δg6, Δneff,1, ⋯, Δneff,5), where σ_neff,i=10^-5 and σ_g,i=5 nm. For both Chebyshev and Butterworth nominal designs, we simulate two sets of design variables with different lower bound (ai, denoted as g_lb here) and upper bound (b_i, denoted as g_ub here), where (g_lb, g_ub)=(0.29, 0.31) and (0.27, 0.33) μm, respectively. The simulation parameters of cases (A–D) are summarized in Table 1. Note that the lower and upper bounds can be, in principle, set differently for each gap, but here they are set to have the same bounds for convenience. Due to space limit, only figures of case A [Chebyshev type, (g_lb, g_ub)=(0.29, 0.31) μm] will be shown in this paper, but all the related statistics are listed in Tables 1 and 2.

Table 1:

Four simulated cases and the associated optimized gaps.

	Filter type	(g_lb, g_ub) (μm)	g→opt (μm)
Case A	Chebyshev	(0.29, 0.31)	(0.3100; 0.3100; 0.3035; 0.2978; 0.3100; 0.3100)
Case B	Chebyshev	(0.27, 0.33)	(0.3054; 0.3300; 0.2858; 0.3014; 0.3300; 0.3300)
Case C	Butterworth	(0.29, 0.31)	(0.2900; 0.3100; 0.2947; 0.2931; 0.3100; 0.2900)
Case D	Butterworth	(0.27, 0.33)	(0.2700; 0.2893; 0.3132; 0.2700; 0.3114; 0.2700)

Table 2:

MSEs of unoptimized nominal design and optimized nominal design under process variations in cases A to D. Each trial has 5000 independent Monte Carlo samples.

	Trial 1			Trial 2			Trial 3
	MSE (nominal)	MSE (optimized)	I (%)	MSE (nominal)	MSE (optimized)	I (%)	MSE (nominal)	MSE (optimized)	I (%)
Case A	0.447	0.402	11.22	0.448	0.403	11.07	0.454	0.406	11.62
Case B	0.447	0.334	33.98	0.448	0.337	32.80	0.454	0.336	35.18
Case C	0.625	0.530	17.83	0.616	0.527	16.83	0.628	0.530	18.46
Case D	0.625	0.512	21.12	0.616	0.509	20.96	0.628	0.502	25.13

The first step is to generate two batches of 5000 Monte Carlo samples with the same fabrication uncertainties, denoted as “training samples” and “test samples”, respectively. Usually 1000 to 5000 samples is a good choice of the training data size. The “training samples” are then used to construct a sparse combined gPC model with total order p=2 using the spgl1 solver [30], [35] in the second step, and the “test samples” are used to verify the capability of the constructed gPC model. The simulated bandwidth pdfs of Monte Carlo samples and gPC model are plotted in Figure 3, showing that the combined gPC model is a good surrogate. Figure 4 plots the bandwidth pdfs of the combined gPC model with fixed design variables g→0, which are randomly selected for the purpose to test the capability of the model. It shows that the combined gPC model can indeed capture accurate statistical information for different design variables and hence can be used as a surrogate in the optimization process.

Figure 3:

3 dB bandwidth pdfs of the combined gPC model with (A) Monte Carlo training samples and (B) Monte Carlo test samples.

$Figure 4: 3 dB bandwidth pdfs of Monte Carlo samples and combined gPC model at given g→0${\vec g_0}$ in case A: (A) g→0=(0.2915, 0.3046, 0.3096, 0.3020, 0.3050, 0.3028)${\vec g_0} = (0.2915,{\rm{ }}0.3046,{\rm{ }}0.3096,{\rm{ }}0.3020,{\rm{ }}0.3050,{\rm{ }}0.3028)$ μm and (B) g→0=(0.295, 0.305, 0.3, 0.305, 0.3, 0.295)${\vec g_0} = (0.295,{\rm{ }}0.305,{\rm{ }}0.3,{\rm{ }}0.305,{\rm{ }}0.3,{\rm{ }}0.295)$ μm.$

Figure 4:

3 dB bandwidth pdfs of Monte Carlo samples and combined gPC model at given g→0 in case A: (A) g→0=(0.2915, 0.3046, 0.3096, 0.3020, 0.3050, 0.3028) μm and (B) g→0=(0.295, 0.305, 0.3, 0.305, 0.3, 0.295) μm.

Our goal is to find the best design whose 3 dB bandwidth is the most robust to (given) process variations. In other words, we try to minimize the MSE of the bandwidth with respect to the nominal bandwidth value (BW₀= 25.6 GHz), and the analytic form of MSE is in Eq. (8). Therefore, we compute MSE in step 3 and then use both the global polynomial optimizer (gloptiPoly3 [29]) and the Matlab global optimization toolbox to solve for global optimum in step 4, as illustrated in Figure 1. The converged solutions for cases A to D are the same with the two solvers and are summarized in Table 1. For all the cases, the construction time of a gPC combined model with design optimization is <15 min (using gloptiPoly3) and 30 min (using Matlab global optimization toolbox).

To verify that a better nominal design is indeed achieved (i.e. its bandwidth is more robust to process variations), we simulate the 3 dB bandwidth of the circuit under process variations with new 5000 Monte Carlo samples for both optimized gaps g→opt and nominal gaps g→nom. The test is performed independently three times, and the corresponding MSEs are summarized in Table 2. It is clearly seen that the circuit with g→opt has smaller MSEs among all the simulated cases. Column I denotes the relative improvement for the optimized design, defined as I=[MSE(nom)-MSE(opt)]/MSE(opt). For case A, the improvement is about 11% in average, whereas case B achieves a larger improvement of about 32.8%–35.2%. Cases A and B are the same filter design (Chebyshev filter), but case B has a larger gap range. This indicates that a better optimum is found in the range considered in case B. The improvement of case C is between 16.8% and 18.5% (larger than Chebyshev filter case A), whereas the MSE improvement in case D (larger gap range) is slightly higher than case C (between 21% and 25.1%) and is smaller compared to the Chebyshev filter (case B).

To visualize the improvement of MSEs, Figure 5 plots the bandwidth pdfs of unoptimized nominal design and optimized design in case A with Trial 1 MC samples. It is shown that the average bandwidth of the unoptimized nominal design is equal to BW₀=25.6 GHz, whereas the average bandwidth of the optimized case is about 25.45 GHz. In addition, the pdf of the optimized design is less dispersed around BW₀. As a result, the MSE improves from about 0.44 (unoptimized nominal design) to 0.40 (optimized design). Note that, although Figures 3–5 all plot the pdfs of bandwidth, they have different meanings. Figure 3 plots the bandwidth pdf where (g→0, Δ→g, n→eff,0) are all varying, whereas Figure 4 refers to designs under fabrication variations with g→0 specified in the caption. Figure 5 also plots the designs under fabrication variations, but g→0 are nominal gaps (g→0=(0.3, 0.3, 0.3, 0.3, 0.3, 0.3) μm) and optimized gaps (g→0=g→opt in case A). Thus, it is natural that their pdfs are all different.

$Figure 5: 3 dB bandwidth pdfs of the unoptimized nominal design (nominal gaps g→0=(0.3, 0.3, 0.3, 0.3, 0.3, 0.3)${\vec g_0} = (0.3,{\rm{ }}0.3,{\rm{ }}0.3,{\rm{ }}0.3,{\rm{ }}0.3,{\rm{ }}0.3)$ μm) and optimized nominal design (optimized gaps g→0=g→opt${\vec g_0} = {\vec g_{opt}}$) with process variations in case A.$

Figure 5:

3 dB bandwidth pdfs of the unoptimized nominal design (nominal gaps g→0=(0.3, 0.3, 0.3, 0.3, 0.3, 0.3) μm) and optimized nominal design (optimized gaps g→0=g→opt) with process variations in case A.

It is interesting to notice that improved bandwidth robustness to fabrication uncertainties comes at the expense of a larger reduction of the average in-band isolation of the filter. Figure 6 plots through-port and drop-port transfer functions of the best (optimized) design in red solid line and the original (unoptimized) nominal design in blue dashed line of case A. It is seen that, for the optimized design, in-band isolation is reduced about 10 dB. This behavior is expected because no constraints were placed on the in-band isolation in the optimization problem (9); hence, its value in the optimized design is not under control. For this reason, there exists a trade-off between improvement in the MSE of bandwidth and reduction of the in-band isolation with respect to the nominal filter, which represents the optimum solution in the ideal situation (without uncertainties), guaranteeing simultaneously the required bandwidth and the best isolation. In our tests, we observed that in case B the isolation of optimized design is about 20 dB lower than nominal design, and the same phenomenon is observed also for cases C and D. This could be improved, for example, requiring in problem (9) the optimized solution to ensure also a minimum value of the isolation as an additional constraint.

Figure 6:

Transfer functions of unoptimized nominal design and optimized design at the through ports and drop ports in case A.

4.3 Further discussion on MSEs and yield

In Section 4.2, it has been demonstrated that the optimized designs are more robust compared to the unoptimized nominal design under fabrication process variations, where σ_neff,i=10^-5 and σ_g,i=5 nm for all four cases. In principle, similar results will be obtained (i.e. optimized designs will outperform the unoptimized nominal design) for different fabrication variations. However, it is interesting to know the capability of the optimized solutions (solved under σ_neff,i=10^-5 and σ_g,i=5 nm) when the fabrication variations change (σ_neff,i≠10^-5 and σ_g,i≠5 nm). This situation is common during real circuit fabrications because the exact uncertainties of the fabrication processes are difficult to predict and may change with time. Of course, if the statistics of fabrication variations are precisely known, then we can always solve (9) to get the optimized design under that fabrication variations as done in Section 4.2.

To perform this test, we calculate the MSE of the nominal design and the case A optimized design under process variations by Monte Carlo simulations, where σ_neff,i is varied from 0.5×10^-5 to 1.5×10^-5 and σ_g,i varied from 1 to 10 nm. The MSE improvement I is reported in Figure 7A. The blue cross in the figure marks the number σ_neff,i=10^-5, σ_g,i=5 nm used in case A, which has I11%. It can be seen that the curves of σ_neff,i=0.5×10^-5 and σ_neff,i=1.5×10^-5 are close to the design case A (σ_neff,i=10^-5), meaning that a change of the effective index variation from 0.5×10^-5 to 1.5×10^-5 has a very small impact on the performance in this example. This is expected because the optical length of the rings has a minor impact on the filter bandwidth at least in the ideal case [33].

Figure 7:

(A) MSE improvement for the optimized design of case A as a function of the σ_g,i and σ_neff,i, where the process variations on the gaps and effective phase indices are assumed to be Gaussian distributed with zero mean and standard deviation σ_g,i and σ_neff,i, respectively. The blue cross represents the MSE improvement for σ_g,i and σ_neff,i used in case A. (B) Yield increment of the optimized designs in cases A–D as a function of M, where 2M is the width of the bandwidth acceptance interval [BW₀-M, BW₀+M].

A better performance of the optimized solution of case A is also observed when gap variation increases (σ_g,i>5 nm) even if the optimized solution here is obtained under the condition σ_g,i=5 nm. The MSE improvement is almost constant at about 11% to 12%. On the contrary, the MSE improvement quickly drops when the gap variation becomes smaller and even negative when σ_g,i<2 nm. This can be explained by the fact that the optimized solution of case A has an average bandwidth slightly different from BW₀ (90 MHz). When the dispersion of the bandwidth pdf decreases due to the reduction of gap variations, its mean value becomes predominant in determining the MSE (i.e. 𝔼[(BW-BW₀)²]). Because the mean bandwidth for the unoptimized nominal design is equal to BW₀, its MSE is almost zero for very small values of σg,i. On the contrary, when σ_g,i increases, the MSE is mainly determined by the dispersion of the 3 dB bandwidth, which is smaller for the optimized solution. At about σg,i=2 nm, the two effects compensate and the nominal and optimized designs have the same performance with respect to the MSE. The above simulation shows that the optimized solution of case A (optimized under gap variation σ_g,i=5 nm) is still effective with gap variation σ_g,i>3 nm; however, a new optimized solution should be computed for smaller gap variation (σg,i<3 nm) to effectively reduce MSEs.

In addition to the MSEs, the yield of a circuit is also a critical index that we would like to observe under process variations. We define a bandwidth acceptance interval [BW₀-M, BW₀+M] as the maximum deviation from the nominal bandwidth (BW₀=25.6 GHz), and the yield is calculated as the integral of the bandwidth pdf in this interval. Because the whole integral of a pdf should be 1, it is obvious that the yield is a number between 0 (i.e. M=0) and 1 (i.e. when M goes to infinite). The higher the yield is, the larger is the number of fabricated filters that meet the specifications.

Figure 7B plots the yield increment (in percentage) for the optimized design compared to the unoptimized nominal design under process variations as a function of the bandwidth acceptance M. For all acceptance intervals, optimized designs outperform unoptimized nominal designs in all four cases. When M is near 0 or larger than 2 GHz, the increment is very small because the yield for both optimized and unoptimized designs approaches 0 and 100%, respectively. On the contrary, when M is between 0.5 and 1 GHz, the maximum increment is observed. The maximum yield increment is about 3.1% (case A, MSE improvement 11%, solid blue line) to 7.3% (case B, MSE improvement 35%, dashed red line). As a clearer visualization of case B, the inset of Figure 7B also plots the yield with respect to M from 0.5 to 1.5 GHz. Case C (dotted black line) has a maximum yield increment similar to case A, whereas case D (dot-dashed magenta line) is about 5.2%.

5 Conclusions

In this paper, we have proposed to construct a sparse combined gPC model for a target performance parameter in photonic circuits. Our sparse model can help solve efficiently a design optimization involving such performance parameter. Specifically, we have applied our methodology to a five-ring coupled resonator filter and designed the gaps so that the 3 dB bandwidth of the optimized designs is more robust to fabrication process variations. The bandwidth MSEs are reduced 11%–35% compared to unoptimized designs, proving the effectiveness of the technique. The yields of the optimized designs are also shown to be higher as an additional benefit.

Acknowledgments

This work was supported through the MIT-SkolTech program, the Progetto Roberto Rocca Seed Funds, the AIM Photonics Center and through the NCN-NEEDS program, which is funded by the National Science Foundation, contract 1227020-EEC, and by the Semiconductor Research Corporation. The authors would also like to thank Dr. Zheng Zhang for his valuable suggestions and comments.

References

[1] Smit M, Leijtens X, Ambrosius H, Bente E, van der Tol J, Smalbrugge B, de Vries T, Geluk E-J, Bolk J, van Veldhoven R, Augustin L, Thijs P, D’Agostino D, Rabbani H, Lawniczuk K, Stopinski S, Tahvili S, Corradi A, Kleijn E, Dzibrou D, Felicetti M, Bitincka E, Moskalenko V, Zhao J, Santos R, Gilardi G, Yao W, Williams K, Stabile P, Kuindersma P, Pello J, Bhat S, Jiao Y, Heiss D, Roelkens G, Wale M, Firth P, Soares F, Grote N, Schell M, Debregeas H, Achouche M, Gentner JL, Bakker A, Korthorst T, Gallagher D, Dabbs A, Melloni A, Morichetti F, Melati D, Wonfor A, Penty R, Broeke R, Musk B, Robbins D. An introduction to InP-based generic integration technology. Semiconduct Sci Technol 2014;29:083001.10.1088/0268-1242/29/8/083001Suche in Google Scholar

[2] Hochberg M, Baehr-Jones T. Towards fabless silicon photonics. Nat Photon 2010;4:492–4.10.1038/nphoton.2010.172Suche in Google Scholar

[3] Chrostowski L, Flueckiger J, Lin C, Hochberg M, Pond J, Klein J, Ferguson J, Cone C. Design methodologies for silicon photonic integrated circuits. In: Proc. SPIE 8989, Smart Photonic and Optoelectronic Integrated Circuits XVI, 89890G, 2014.10.1117/12.2047359Suche in Google Scholar

[4] Melati D, Morichetti F, Canciamilla A, Roncelli D, Soares FM, Bakker A, Melloni A. Validation of the building-block-based approach for the design of photonic integrated circuits. J Lightwave Technol 2012;30:3610–6.10.1109/JLT.2012.2223658Suche in Google Scholar

[5] Chen X, Mohamed M, Li Z, Shang L, Mickelson AR. Process variation in silicon photonic devices. Appl Opt 2013;52:7638–47.10.1364/AO.52.007638Suche in Google Scholar PubMed

[6] Melati D, Alippi A, Melloni A. Waveguide-based technique for wafer-level measurement of phase and group effective refractive indices. J Lightwave Technol 2015;34:1293–9.10.1109/JLT.2015.2500919Suche in Google Scholar

[7] Selvaraja SK, Bogaerts W, Dumon P, Van Thourhout D, Baets R. Subnanometer linewidth uniformity in silicon nanophotonic waveguide devices using CMOS fabrication technology. Select Top Quantum Electron IEEE J 2010;16:316–24.10.1109/JSTQE.2009.2026550Suche in Google Scholar

[8] Liew SF, Ge L, Redding B, Solomon GS, Cao H. Pump-controlled modal interactions in microdisk lasers. Phys Rev A 2015;91:043828.10.1103/PhysRevA.91.043828Suche in Google Scholar

[9] Melati D, Melloni A, Morichetti F. Real photonic waveguides: guiding light through imperfections. Adv Opt Photon 2014;6:156–224.10.1364/AOP.6.000156Suche in Google Scholar

[10] Cheung S, Su OT, Yoo K. Ultra-compact silicon photonic 512 × 512 25 GHz arrayed waveguide grating router. Select Top Quantum Electron IEEE J 2014;20:310–6.10.1109/JSTQE.2013.2295879Suche in Google Scholar

[11] Chrostowski L, Wang X, Flueckiger J, Wu Y, Wang Y, Fard ST. Impact of fabrication non-uniformity on chip-scale silicon photonic integrated circuits. In: Optical Fiber Communication Conference, OSA Technical Digest (online) (Optical Society of America, 2014), paper Th2A.37.10.1364/OFC.2014.Th2A.37Suche in Google Scholar

[12] Melati D, Lovati E, Melloni A. Statistical process design kits: analysis of fabrication tolerances in integrated photonic circuits. In: Integrated photonics Research, Silicon and Nanophotonics. Optical Society of America, 2015:IT4A–5.10.1364/IPRSN.2015.IT4A.5Suche in Google Scholar

[13] Wang X, Shi W, Yun H, Grist S, Jaeger NAF, Chrostowski L. Narrow-band waveguide Bragg gratings on SOI wafers with CMOS-compatible fabrication process. Opt Express 2012;20:15547–58.10.1364/OE.20.015547Suche in Google Scholar

[14] Bogaerts W, Fiers M, Dumon P. Design challenges in silicon photonics. Sel Top Quantum Electron IEEE J 2014:20:1–8.10.1109/JSTQE.2013.2295882Suche in Google Scholar

[15] Cassano D, Morichetti F, Melloni A. Statistical analysis of photonic integrated circuits via polynomial-chaos expansion. In: Advanced Photonics 2013, OSA Technical Digest (online) (Optical Society of America, 2013), paper JT3A.8.10.1364/IPRSN.2013.JT3A.8Suche in Google Scholar

[16] Ghanem R, Spanos PD. A stochastic Galerkin expansion for nonlinear random vibration analysis. Probabilist Eng Mech 1993;8:255–64.10.1016/0266-8920(93)90019-RSuche in Google Scholar

[17] Zhang Z, El-Moselhy TA, Elfadel IM, Daniel L. Stochastic testing method for transistor-level uncertainty quantification based on generalized polynomial chaos. Comput Aided Des Integr Circuits Syst IEEE Trans 2013;32:1533–45.10.1109/TCAD.2013.2263039Suche in Google Scholar

[18] Xiu D, Karniadakis GE. Modeling uncertainty in flow simulations via generalized polynomial chaos. J Comput Phys 2003;187:137–67.10.21236/ADA461813Suche in Google Scholar

[19] Candès EJ, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. Inf Theory IEEE Trans 2006;52:489–509.10.1109/TIT.2005.862083Suche in Google Scholar

[20] Chen SS, Donoho DL, Saunders MA. Atomic decomposition by basis pursuit. SIAM J Sci Comput 1998;20:33–61.10.1137/S1064827596304010Suche in Google Scholar

[21] Doostan A, Owhadi H. A non-adapted sparse approximation of pdes with stochastic inputs. J Comput Phys 2011;230:3015–34.10.1016/j.jcp.2011.01.002Suche in Google Scholar

[22] Zhang Z, Weng T-W, Daniel L. A big-data approach to handle process variations: uncertainty quantification by tensor recovery. arXiv:1603.06119 (2016) [arXiv preprint].10.1109/SaPIW.2016.7496314Suche in Google Scholar

[23] Dodson M, Parks GT. Robust aerodynamic design optimization using polynomial chaos. J Aircraft 46:635–46, 2009.10.2514/1.39419Suche in Google Scholar

[24] Zhao L, Dawes WN, Parks G, Jarrett JP, Yang S. Robust airfoil design with respect to boundary layer transition. In: Proceedings of 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, vol. 4, 2009.10.2514/6.2009-2273Suche in Google Scholar

[25] Ghisu T, Jarrett JP, Parks GT. Robust design optimization of airfoils with respect to ice accretion. J Aircraft 2011;48:287–304.10.2514/1.C031100Suche in Google Scholar

[26] Ghisu T, Parks GT, Jarrett JP, Clarkson PJ. Robust design optimization of gas turbine compression systems. J Propul Power 2011;27:282–95.10.2514/1.48965Suche in Google Scholar

[27] Eldred MS. Design under uncertainty employing stochastic expansion methods. Int J Uncertain Quantif 201;1:119–46.10.2514/6.2008-6001Suche in Google Scholar

[28] Adams BM, Ebeida MS, Eldred MS, Jakeman JD, Swiler LP, Stephens JA, Vigil DM, Wildey TM, Bohnhoff WJ, Eddy JP, Hu KT, Bauman LE, Hough PD. Dakota, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis. Technical report, Sandia National Laboratories (SNL-NM), Albuquerque, NM (United States), 2014.10.2172/1177077Suche in Google Scholar

[29] Henrion D, Lasserre J-B, Löfberg J. Gloptipoly 3: moments, optimization and semidefinite programming. Optim Method Softw 2009;24:761–79.10.1080/10556780802699201Suche in Google Scholar

[30] van den Berg E, Friedlander MP. SPGL1: a solver for large-scale sparse reconstruction, 2007. http://www.cs.ubc.ca/labs/scl/spgl1.Suche in Google Scholar

[31] Lasserre JB. A semidefinite programming approach to the generalized problem of moments. Math Program 2008;112:65–92.10.1007/s10107-006-0085-1Suche in Google Scholar

[32] Melloni A, Martinelli M. Synthesis of direct-coupled-resonators bandpass filters for wdm systems. J Lightwave Technol 2002;20:296–303.10.1109/50.983244Suche in Google Scholar

[33] Madsen CK, Zhao JH. Optical filter design and analysis: a signal processing approach. 1st ed. New York, NY, USA, John Wiley & Sons, Inc., 1999.10.1002/0471213756Suche in Google Scholar

[34] Aspic. http://aspicdesign.com. Accessed 03-01-2016.Suche in Google Scholar

[35] van den Berg E, Friedlander MP. Probing the Pareto frontier for basis pursuit solutions. SIAM J Sci Comput 2008;31: 890–912.10.1137/080714488Suche in Google Scholar

Received: 2016-3-12

Revised: 2016-5-8

Accepted: 2016-5-21

Published Online: 2016-7-30

Published in Print: 2017-1-6

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.

Artikel in diesem Heft

https://doi.org/10.1515/nanoph-2016-0110

Schlagwörter für diesen Artikel

design optimization; generalized polynomial chaos; process variations; photonic integrated circuits; yield; coupled ring resonator filters

Creative Commons

BY-NC-ND 3.0