1. Introduction

As silicon photonic integration technology rapidly advances, a growing body of research in silicon photonics is moving from academic laboratories to the industry, where manufacturing yields become a critical concern. Due to the high refractive index contrast, the performance of silicon-on-insulator-based photonic integrated circuits (PICs) is more vulnerable to manufacturing process variations compared to PICs on alternative material platforms [1–3]. Hence, it is imperative to undertake a quantitative analysis of how process variations might impact the performance of PICs prior to manufacturing. This analysis is commonly referred to as variability analysis or uncertainty quantification.

The manufacturing process has a direct impact on the geometric and physical design parameters in PICs. For example, the width of waveguides may vary between dies and wafers due to factors like exposure and etching processes, as well as the local pattern density of the mask layout. Similarly, slight alterations in ion implantation conditions, such as the implanted dose, ion beam angle, and energy, can introduce non-uniform and stochastic doping concentrations. The stochastic nature of these parameters can result in notable deviations in the optical properties of photonic devices, such as effective index and coupling coefficient. When these devices are interconnected to form a circuit, the cumulative impact of fabrication variations can ultimately degrade the overall performance of the circuit.

The conventional method for assessing the performance variations of PICs is the Monte Carlo (MC) method. This involves running numerous simulations to gather statistical information. Consequently, it could impose a substantial computational burden when a single simulation is computationally expensive. An alternative solution is to leverage the Polynomial Chaos (PC) expansion technique, which is able to represent stochastic quantities as the sum of weighted orthogonal polynomial basis functions [4]. Many techniques based on the PC method have been recently developed to quantify the variations in frequency response of passive PICs [5–8]. Considering the time-domain evaluation of PICs is more intuitive and informative, the author has previously suggested a time-domain stochastic modeling approach for passive PICs [9]. This approach allows for the accurate assessment of time-domain statistical information of passive PICs via simulating the built model with any given input signals. While this method is accurate, it comes with certain limitations. In the modeling process, deterministic augmented scattering matrices are formulated through suitable combinations of the PC coefficients derived from the stochastic scattering matrices of the original PICs. The size of these augmented matrices can become considerably large, particularly when a large number of basis functions is utilized. Consequently, building models for the large-size augmented scattering matrices while guaranteeing accuracy and passivity could be exceedingly challenging. Even if successfully managed, it results in significantly large time-domain state-space models, leading to time-consuming circuit simulations.

To address this challenge, we develop an easy-to-implement and efficient time-domain variability analysis method for passive PICs. The method leverages PC expansion technique and builds deterministic models for stochastic PICs. The models are constructed by properly interconnecting the passive PICs defined at a selected set of random parameter samples and their port signals are the PC coefficients of the stochastic port signals of the studied PICs. Considering that there is currently no standard modeling language for PICs, the constructed models are implemented using Verilog-A, the de facto standard language for defining and distributing models in the semiconductor industry. This choice enables the models to be simulated across various sophisticated commercial software platforms and co-simulated with electronic integrated circuits [10,11]. Finally, with a given input, the time-domain statistical information of the PIC under study can be computed through a single simulation of these Verilog-A models.

The structure of this paper is as follows: Section 2 provides a brief introduction about PC expansion of the signal and system of passive PICs. In Section 3, we review the existing stochastic modeling method and modify it to improve efficiency. Section 4 proposes the new stochastic modeling method and demonstrates how to implement the built model in Verilog-A. The proposed method is validated in Section 5 through the statistical simulations of two wavelength filters. Finally, conclusions are drawn in Section 6.

2. Polynomial chaos expansion of passive PICs

The frequency-dependent behaviors of passive PICs are fully described through their scattering (S) matrices, which can be acquired via measurements or electromagnetic simulations in the optical frequency range of interest. In the presence of random variations, the relationship between the port signals and S matrices is

According to the PC theory, the stochastic quantities ${\boldsymbol {a}}({\boldsymbol {\xi }},f)$, ${\boldsymbol {b}}({\boldsymbol {\xi }},f)$, and ${\boldsymbol {S}}({\boldsymbol {\xi }},f)$ can be expanded as the sum of a series of weighted orthonormal basis functions

3. Review and improvement of the existing stochastic modeling method for passive PICs

3.1 Review of the existing time-domain stochastic modeling method

This section reviews the existing time-domain stochastic modeling method for passive PICs that was proposed in the previous work [9]. The technique starts with substituting the stochastic quantities in (1) with their PC expansions in (2)

Then, projecting (3) onto the $p-$th basis function via Galerkin projection gives

Now, (5) describes a new deterministic system represented by ${\boldsymbol {S}}_{PC}(f)$, whose port signals are the PC coefficients of the port signals of the original stochastic system: ${\boldsymbol {S}}_{PC}(f)$ is $M+1$ times larger than the original system under study in terms of ports number. At this point, a time-domain state-space model can be computed from ${\boldsymbol {S}}_{PC}(f)$ via the Vector Fitting (VF) technique [12] and then shifted to baseband by carrier frequency for efficient time-domain simulations [13,14]. The ‘shifting’ operation is to derive a baseband representation for the PICs, enabling the adoption of relatively large time steps in time-domain simulations.

3.2 Improvement of the existing time-domain stochastic modeling method

The aforementioned method can be improved by incorporating the new baseband modeling technique proposed in [15], known as Complex Vector Fitting (CVF) method. With this method, the S matrices ${\boldsymbol {S}}({\boldsymbol {\xi }},f)$ defined in optical frequency range are first shifted to baseband by the carrier frequency, leading to baseband S matrices ${\boldsymbol {S}}({\boldsymbol {\xi }},\bar {f})$. Note that, baseband S matrices are complex envelope representation of the original system, which are non-physical and do not comply with the Hermitian symmetric constraint. Meanwhile, the forward and backward optical waves can also be represented with their baseband equivalents. For more details about the definitions of baseband equivalent signals and the system of passive PICs, please refer to [15].

Then, the same procedure, as detailed from (3) to (6), is applied to ${\boldsymbol {S}}({\boldsymbol {\xi }},\bar {f})$, yielding the corresponding large deterministic augmented system denoted as ${\boldsymbol {S}}_{PC}(\bar {f})$. Similar to ${\boldsymbol {S}}({\boldsymbol {\xi }},\bar {f}),$ ${\boldsymbol {S}}_{PC}(\bar {f})$ does not exhibit Hermitian symmetry. As a result, the VF technique is not applicable to modeling ${\boldsymbol {S}}_{PC}(\bar {f})$ any more. Instead, the CVF method can be adopted in this scenario to construct a time-domain state-space model for ${\boldsymbol {S}}_{PC}(\bar {f})$

The primary distinction between the improved and the existing methods is as follows: In [9], S matrices are initially modeled in the optical range, and subsequently, the corresponding time-domain model constructed through the VF technique is shifted to the baseband for efficient time-domain simulations. In contrast, in this section, the S matrices are first shifted to the baseband and then directly modeled via the CVF method. The key advantage is that the enhanced method can generate a state-space model with half the size for the same applications, facilitating more efficient time-domain simulations. Furthermore, the substantial size of the augmented S matrices in this context often leads to significantly large state-space models. Thus, the efficiency improvement is more pronounced as the size of these large models is reduced by half. Interested readers are encouraged to refer to [9,15] for more details.

Once the model (7) is computed, a single time-domain simulation with any given inputs generates all the PC coefficients for the port signals of the stochastic PICs. Unfortunately, despite the superiority of the improved method over the one presented in [9], there are still several limitations. Firstly, calculating both ${\boldsymbol {S}}_{PC}(f)$ and ${\boldsymbol {S}}_{PC}(\bar {f})$ requires to compute $\langle \varphi _{j}({\boldsymbol {\xi }})\varphi _{k}({\boldsymbol {\xi }}),\varphi _{p}({\boldsymbol {\xi }})\rangle$, which is analytically available for standard polynomial classes, but is non-trivial to calculate for arbitrary or non-Gaussian correlated distributions. Secondly, modeling ${\boldsymbol {S}}_{PC}(\bar {f})$ using the CVF method while guaranteeing accuracy and passivity at the same time could be exceedingly challenging when the number $M+1$ of basis functions is large. Even if successfully managed, it leads to a significantly large state-space model (7) [9]. Implementing and simulating such a large and indivisible model in Verilog-A poses substantial challenges. To overcome these issues, we propose a new stochastic modeling approach in the next section.

4. Proposed Verilog-A based variability analysis method for passive PICs

4.1 New stochastic modeling method of passive PICs

Assuming the baseband S matrices of the passive PIC under study is ${\boldsymbol {S}}({\boldsymbol {\xi }},\bar {f})$, the CVF method can construct a time-domain state-space model from ${\boldsymbol {S}}({\boldsymbol {\xi }},\bar {f})$ in the form

According to the work [16,17], the $M+1$ PC coefficients of the stochastic quantities ${\boldsymbol {a}}({\boldsymbol {\xi }},t)$ can be accurately determined with the values of ${\boldsymbol {a}}({\boldsymbol {\xi }},t)$ at $M+1$ selected samples $\{{\boldsymbol {\xi }}_k\}_{k=0}^{M}$ for independent random variables. For non-Gaussian correlated random variables, the number of samples required is larger than $M+1$ and is denoted as $K+1$ [18]. In the rest of the paper, $K \geqslant M$ is used to encompass different cases. At each sample ${\boldsymbol {\xi }}_k$, the S matrices ${\boldsymbol {S}}({\boldsymbol {\xi }}_k,f)$ are evaluated in the frequency range of interest and the state-space matrices ${\boldsymbol {a}}({\boldsymbol {\xi }}_k)$, ${\boldsymbol {b}}({\boldsymbol {\xi }}_k)$, ${\boldsymbol {C}}({\boldsymbol {\xi }}_k)$, and ${\boldsymbol {D}}({\boldsymbol {\xi }}_k)$ can be obtained via the CVF modeling method. By representing the forward and backward waves with their PC expansion, the CVF model at sample ${\boldsymbol {\xi }}_k$ becomes

Considering that there are $K+1$ models for the $K+1$ samples, we formulate (10) at $K+1$ samples into matrix form

 

Fig. 1. The structure of the proposed model for time-domain variability analysis of passive PICs.

Download Full Size | PDF

When multiple devices and subcircuits in a circuit experience variations, their models in the form shown in Fig. 1 can be hierarchically connected in accordance with the original network topology. The PC coefficients of the signals at all ports and internal nodes can be obtained via a single time-domain simulation. It is worth mentioning that while this method is proposed for passive PICs, it has the potential for extension to general PICs containing both passive and active components. In particular, active components such as modulators can be readily incorporated into the framework depicted in Fig. 1 by substituting the CVF models with the models of the active components.

4.2 Verilog-A realization of the proposed model

This section demonstrates how to implement the model (11) in Verilog-A. The signals and models presented in this work are baseband representations, thus they are complex-valued. However, Verilog-A does not inherently support complex numbers. To solve this issue, a two-signal bus is used to represent the real and the imaginary parts of baseband signals, as has been done in [11]. The CVF models contained in Fig. 1 are also complex-valued, but they can be converted into real-valued form by separating the real and imaginary parts of the complex signals and matrices in (8) according to [15].

The model (11) primarily involves three mathematical operations: addition, multiplication, and solving the real-valued state-space models derived from the CVF models. Note that, the matrix inverse can be precomputed in any other software, such as Matlab. Therefore, the main focus is on how to implement the real-valued state-space model in Verilog-A. Considering that the state-space models are inherently first-order ordinary differential equations, the derivative operation in Verilog-A is well suited to be utilized. For example, a CVF model with one input, one output, and two state variables is expressed as

 

Fig. 2. Verilog-A realization of the real-valued state-space model (15).

Download Full Size | PDF

The discipline ‘optical’ in Fig. 2 is defined in a manner consistent with the presentation in [11]. Thanks to the high-level abstraction provided by Verilog-A, implementing the model (15) is both concise and efficient. Moreover, the transformation of the model into Verilog-A code can be readily automated using any programming language, with Matlab being the tool of choice in our case. Once the Verilog-A realization is generated, it can be simulated and reused across various commercial software platforms that support Verilog-A.

5. Validation and numerical results

This section validates the proposed modeling and simulation technique with two application examples. The models are computed and then converted to Verilog-A code in Matlab while the Verilog-A code of the models are simulated in Cadence Virtuoso on a personal computer with Intel Core i7 processor and 32 GB RAM.

5.1 Variability analysis of a three-ring resonator filter

The section focuses on the study of a three-ring resonator filter, illustrated in Fig. 3(a). The rings are configured in a race-track shape, with ring radius specified as 10 $\mu$m and waveguide width set at 450 nm and a loss of 2 dB/cm. The gaps ($g_1$ and $g_4$) between the access waveguide and the ring are designed to be 200 nm, while inter-ring gaps ($g_2$ and $g_3$) are set at 120 nm. To account for fabrication variations, it is assumed that all gaps vary independently and follow a Gaussian distribution with a standard deviation of 2 nm. The objective is to quantify the impact of these variations on the time-domain performance of the filter given a specific input sequence. It’s important to note that a slightly low standard deviation of the gap is used here to prevent the generation of a larger model by the improved version of the existing method, which would pose significant challenges for simulation in Cadence Virtuoso.

 

Fig. 3. (a) The structure of the three-ring resonator filter, (b) Variations of frequency responses of the three-ring resonator filter.

Download Full Size | PDF

The filter is designed and simulated in IPKISS to obtain the S matrices at a set of frequency points in the range [192.17, 194.17] THz. Due to variations in the gaps $g_1$ to $g_4$, the frequency responses at the pass and drop port also exhibit variations, as illustrated in Fig. 3(b). The extensive collection of light blue and red lines represents the frequency responses obtained from 10,000 MC samples of $\{{\boldsymbol {g}}_k=[g_1, g_2, g_3, g_4]_k\}_{k=1}^{10000}$, whereas the blue and red lines correspond to the nominal case. As the S parameters have been shifted to baseband relative to the center frequency of the passband/stopband, the frequency axis in Fig. 3(b) is centered around zero.

The time-domain variability analysis is conducted using three methods: the MC method, the improved version of the existing method (referred to as the improved method), and the proposed method. In the MC method, 10,000 CVF models with 8 poles are built from the S matrices ${\boldsymbol {S}}({\boldsymbol {g}}_k)$ of the filter at 10,000 of ${\boldsymbol {g}}_k$. For the improved method, the PC order $P$ is selected to be two for the four random variables, resulting in a total of 15 PC coefficients. The ${\boldsymbol {S}}_{PC}$ at different frequency points are calculated via linear regression with ${\boldsymbol {S}}({\boldsymbol {g}}_k)$ at $\{{\boldsymbol {g}}_k\}_{k=1}^{15}$. Notably, the 15 samples of ${\boldsymbol {g}}_k$ are chosen using the method outlined in [9,16]. Subsequently, the CVF modeling method is employed to build the state-space model (7) for ${\boldsymbol {S}}_{PC}$ with 12 poles. It is worth mentioning that the size of ${\boldsymbol {S}}_{PC}$ is 60 by 60, while the state-space matrices ${\boldsymbol {a}}_{PC}\in \mathbb {C}^{720\times 720}$, ${\boldsymbol {b}}_{PC}\in \mathbb {R}^{720\times 60}$, ${\boldsymbol {C}}_{PC}\in \mathbb {C}^{60\times 720}$, and ${\boldsymbol {D}}_{PC}\in \mathbb {C}^{60\times 60}$. Finally, the proposed model in the form depicted in Fig. 1 is also constructed with the same 15 samples. All the models are then converted into Verilog-A realizations prepared for simulation in Cadence Virtuoso.

The time-domain performance is analyzed when exciting the filter at the input port with a modulated 50 Gbit/s PRBS11 (Pseudo-Random Binary Sequence 11) signal featuring a smooth waveform and 2047 symbols. The time-domain simulations are then carried out with the three methods in Cadence Virtuoso. The mean and standard variations of the signals at drop and pass ports can be immediately computed from the outputs of the improved method and the proposed method. Specifically, ${{\boldsymbol {b}}_0}(t)$ represents the mean of ${\boldsymbol {b}}({\boldsymbol {\xi }},t)$, while the variance can be estimated with the remaining coefficients through ${\sum _{m=1}^{M}{{\boldsymbol {b}}_m^2}(t)}$. Meanwhile, the mean and variance of the MC method can also be calculated by analyzing the 10,000 simulation results, serving as a benchmark. Figure 4 compares the mean and standard deviation (sigma) of the outputs at the drop port from different methods. Since the filter is modeled as baseband equivalent systems, the outputs are complex numbers even if the input is real. Therefore, the real and imaginary parts are plotted in the figure. It is observed that the presence of Gaussian random variables in the filter leads to non-Gaussian time-domain outputs, as evidenced by the fact that the $mean\pm 3sigma$ range does not encompass the results from the 10,000 MC simulations. The results in Fig. 4 also demonstrate the improved method and the proposed method are as accurate as the MC method.

 

Fig. 4. Variation of the output at the drop port of the three-ring resonator filter. The gray lines are the results from 10,000 MC simulations. The blue, red and green lines correspond to the $mean+3sigma$, $mean$, and $mean-3sigma$ from the MC method, respectively. The (blue, red, and green) markers $\blacktriangle$ and $\bullet$ denote the same quantities from the improved method and the proposed method, respectively.

Download Full Size | PDF

To facilitate a more detailed comparison between the outputs corresponding to the PRBS11 input from the proposed method and the MC method, eye diagrams of the drop-port signals are displayed in Fig. 5. Note that, with the PC coefficients obtained from the proposed method, the drop-port outputs at 10,000 samples of ${\boldsymbol {g}}$ can be immediately calculated. This visualization further confirms the accuracy of the proposed method. Given that the improved method achieves a comparable level of accuracy to the proposed method, its eye diagram is omitted in Fig. 5. For a detailed characterization of the eye diagrams, we additionally calculate the probability density function (PDF) for the eye height variations caused by the variances in the gaps, as shown in Fig. 6.

 

Fig. 5. Eye diagrams. (a) PRBS11 input sequence, (b) Nominal output at drop port, (c) 10,000 outputs at drop port calculated from the proposed method, (d) 10,000 outputs at drop port from the MC method.

Download Full Size | PDF

 

Fig. 6. PDF of the eye heights calculated from the MC method and the proposed method.

Download Full Size | PDF

Beside accuracy, efficiency is another main concern for variability analysis methods. Therefore, we present the modeling time and simulation time of the three methods in Table 1. It demonstrates that the proposed method is about 11 times faster than the improved method and 800 times faster than the MC method.

Table 1. Efficiency of the proposed technique

View Table | View all tables in this article

5.2 Variability analysis of a MZI-based lattice filter

In this section, a MZI-based lattice filter is designed, as shown in Fig. 7(a), and simulated in IPKISS to obtain its S matrices. All waveguides in the filter have a width of 500 nm and thickness of 220 nm. The coupling gap of the nine directional couplers are set to be 200 nm, resulting in a same per-length coupling coefficient $\kappa$ of 0.0472. The cross-coupled power is tuned by adjusting the coupling length $L$. Due to gap variations introduced by manufacturing process, $\kappa$ can vary in a certain range [19]. Given that the nine directional couplers are situated at different but nearby locations on wafers, we assume that $\kappa _1$ to $\kappa _9$ are non-Gaussian correlated random parameters. According to [18], it is possible to represent the distribution of ${\boldsymbol {\kappa }}=[\kappa _1,\ldots,\kappa _9]^T$ with a Gaussian mixture model

 

Fig. 7. (a) The structure of the MZI-based lattice filter, (b) Variation of transmissions in different channels of the MZI-base lattice filter.

Download Full Size | PDF

To construct the proposed model, a PC order of $P=2$ is selected, leading to a total of $M+1=55$ PC coefficients, taking into account that the number of random variables is $d=9$. Employing the technique proposed in [18], the corresponding orthonormal basis functions are computed, and a subset of 83 samples is carefully selected to ensure the accurate construction of the model. In this case, the number of samples, denoted as $K+1$, exceeds the number of basis functions, $M+1$. Therefore, the Moore-Penrose pseudo-inverse is employed in the formulation given by (11). Subsequently, the proposed model, depicted in the structure shown in Fig. 1, is implemented in Verilog-A for time-domain simulations. Concurrently, the CVF models for the filter at 10,000 MC samples are constructed and converted into Verilog-A representation. In this example, the improved method is not demonstrated due to the difficulties mentioned in Section 3.2. Indeed, calculating ${\boldsymbol {S}}_{PC}$ requires the evaluation of $\langle \varphi _{j}({\boldsymbol {\xi }})\varphi _{k}({\boldsymbol {\xi }}),\varphi _{p}({\boldsymbol {\xi }})\rangle$ which is a non-trivial task when dealing with random variables characterized by non-standard distributions. In addition, the dimensionality of ${\boldsymbol {S}}_{PC}$ will expand to 275 by 275 when restricting the analysis to five ports utilizing 55 basis functions. Consequently, the corresponding CVF model will also be substantial, rendering it extremely difficult, if not unfeasible, to simulate such a large and indivisible model in Cadence Virtuoso.

In order to assess the time-domain performance of the filter, we generate a PAM4 signal characterized by a raised cosine pulse shape with a symbol rate of 80 Gbaud/s and a total length of 2000 symbols. Assuming the modulation of this signal onto the carrier corresponding to channel one and its propagation into port P1, the task is to analyze variations in both the transmission at P5 and crosstalk at P2 to P4. Time-domain simulations of both the proposed model and the 10,000 CVF models are carried out using Cadence Virtuoso. The mean and standard variations of the outputs at P2 to P5 can be analytically computed from the outputs of the proposed model. With the MC method, the same quantities can also be computed using the simulation results of the 10,000 CVF models, serving as a benchmark.

Given that the outputs are complex-valued signals, Fig. 8(a) compares the $mean$ and $mean\pm 3sigma$ of the real and imaginary components of the output at port P5. These values are computed using both the MC method and the proposed method. Due to inadequate out-of-band rejection of the channels two to four, a fraction of the input signal designated for channel one passes through to channels two to four, resulting in what is considered as crosstalk. For instance, Fig. 8(b) illustrates the variation of crosstalk in the channel two. Note that the figures only depict outputs from a small time window for clarity, as the entire sequence signal is quite lengthy. For a more comprehensive representation of the accuracy of the proposed model, we compare the computed mean and standard deviation in Fig. 9(a) and Fig. 9(b), respectively. All of these results demonstrate that the proposed method performs with accuracy comparable to the MC method.

 

Fig. 8. (a) Variation of the transmission at port P5 of the MZI-base lattice filter, (b) Variation of the crosstalk at port P2 of the MZI-base lattice filter. The gray lines are the results from 10,000 MC simulations. The $mean$ (red color), $mean+3sigma$ (blue color), and $mean-3sigma$ (green color) calculated from the MC results (lines) and the proposed method (dots) are compared.

Download Full Size | PDF

 

Fig. 9. (a) Mean of the transmission at port P5 of the MZI-base lattice filter, (b) Standard deviation of the transmission at port P5 of the MZI-base lattice filter.

Download Full Size | PDF

It is intriguing to note, upon closer examination of the zoomed window in Fig. 8(a), that the $mean\pm 3sigma$ curves do not encompass the 10,000 MC results. This suggests that the outputs are not Gaussian distributed, where $mean\pm 3sigma$ would typically include 99.7 percent of the data. To verify this observation, we present the PDF of the output at the time instant $t=3.852$ ns in Fig. 10. The PDF results of the proposed model are derived by calculating 10,000 MC samples using the obtained PC coefficients of the signals at that specific time instant. This reveals that the $mean\pm 3sigma$ range does not align with the distribution of the data, which is consistent with the plot in Fig. 8(a). Once again, the comparison of the PDF plots validates the accuracy of the proposed method.

 

Fig. 10. PDF of the transmission ($t=3.852$ ns) at port P5 of the MZI-base lattice filter. Results from the proposed method is compared to that from the MC method while the $mean\pm 3sigma$ range is also indicated.

Download Full Size | PDF

Figures 8, 9, and 10 focus on comparing results within a limited time window. To extend the comparison across the entire output sequences, we construct eye diagrams for the input, nominal transmission, and the 10,000 transmissions from both the MC and proposed methods, as depicted in Fig. 11. Notably, it is observed that the eye size varies significantly under stochastic effects, and the eye diagram from the proposed method offers a highly accurate estimation of the MC results.

 

Fig. 11. Eye diagrams. (a) PAM4 input sequence at P1, (b) Nominal transmission at P5, (c) 10,000 transmissions from the proposed method, (d) 10,000 transmissions from the MC method.

Download Full Size | PDF

The primary advantages of the proposed method lie in its ability to deliver comparable accuracy to the MC method while exhibiting significantly higher efficiency. The modeling time and simulation time for both methods are detailed in Table 2, illustrating that the proposed method achieves a speedup of over 100 times compared to the MC method throughout the entire modeling and simulation process. It’s worth noting that the first example showcases a much higher speedup compared to this example. The reason for this difference is that the proposed model in the first example, which involves four random variables, is constructed using only 15 samples. In contrast, in this example, which involves nine random variables, 83 samples are employed for model construction.

Table 2. Efficiency of the proposed technique

View Table | View all tables in this article

6. Conclusion

This paper proposed a circuit-level variability analysis method for passive PICs with the potential for extension to general PICs. It is applicable to PICs with random parameters that follow both standard and non-standard distributions, mirroring real-world scenarios. Furthermore, this approach is considerably easy to implement and exhibits superior scalability compared to existing method and the improved method, particularly as the numbers of ports and random parameters in the studied PICs increase. Lastly, the proposed model has been implemented in Verilog-A, facilitating simulations across various commercial simulators in semiconductor industry and enabling co-simulation with electronic integrated circuits. The accuracy and efficiency of the proposed method, compared to the improved method and the MC method, are validated through two application examples.