Open Access
Add to BibSonomyAbstract
This work develops an enhanced Monte Carlo (MC) simulation methodology to predict the impacts of layout-dependent correlated manufacturing variations on the performance of photonics integrated circuits (PICs). First, to enable such performance prediction, we demonstrate a simple method with sub-nanometer accuracy to characterize photonics manufacturing variations, where the width and height for a fabricated waveguide can be extracted from the spectral response of a racetrack resonator. By measuring the spectral responses for a large number of identical resonators spread over a wafer, statistical results for the variations of waveguide width and height can be obtained. Second, we develop models for the layout-dependent enhanced MC simulation. Our models use netlist extraction to transfer physical layouts into circuit simulators. Spatially correlated physical variations across the PICs are simulated on a discrete grid and are mapped to each circuit component, so that the performance for each component can be updated according to its obtained variations, and therefore, circuit simulations take the correlated variations between components into account. The simulation flow and theoretical models for our layout-dependent enhanced MC simulation are detailed in this paper. As examples, several ring-resonator filter circuits are studied using the developed enhanced MC simulation, and statistical results from the simulations can predict both common-mode and differential-mode variations of the circuit performance.
© 2017 Optical Society of America
1. Introduction
The high refractive index contrast of silicon-on-insulator (SOI) photonics enables tight confinement of light in sub-micrometer waveguides and sharp waveguide bends, making SOI promising for developing photonic integrated circuits (PICs) with high integration densities. In applications such as wavelength division multiplexing (WDM), PICs often require precise matching of the central wavelength and the waveguide propagation constant between components on a chip (e.g., ring modulators, optical filters). However, manufacturing variability is a challenge [1] in PIC designs. Fabrication errors in waveguide width and height have significant impacts on the propagation constant of light in waveguides; therefore, circuits such as interferometers are highly sensitive to manufacturing. For instance, Mach-Zehnder interferometers are more sensitive to differential variations on the two phase arms than the common-mode variations to the entire device [2]. For designers, it is critical to quantify fabrication variations and take such correlated physical variations into account in the PIC designs. With estimation on the performance variations of PICs, appropriate strategies (e.g., thermal tuning) can be developed and implemented to compensate such variations, and the cost implications for such compensation strategies (e.g., power consumption) can also be determined.
Characterizing manufacturing variations in a wafer scale is challenging. Although the width and height of fabricated SOI waveguides can be characterized by using scanning electron microscope (SEM) imaging [3] and atomic force microscope (AFM) mapping [4], respectively, they are rather costly and time-consuming. A more efficient method to characterize fabrication variations is investigating the spectral response variations of identical SOI devices (e.g., micro-disk resonators [5] and Bragg gratings [6, 7]), and using analytic algorithms to extract the waveguide width and height from such spectral response variations. However, most demonstrated methods require either dual polarizations measurements [5] or complex spectrum reflection measurements [6]. More efforts are required to simplify the characterization methods.
Taking spatially-dependent manufacturing variations into account in PIC designs is another challenge for designers. In electronic designs, corner analysis and Monte Carlo (MC) analysis are typically used to study the impacts of manufacturing variability. The corner analysis only predicts the performance at the process corners of the each design parameter, without exploring the cases having parameters falling within the range of process corners. Unfortunately, this approach is not appropriate for photonics [2]. The simplest MC approach applies the same random variation to all of the circuit components, and thus, simulations consider common-mode variability only and are not able to capture correlated variations between components. The MC approach can also assign differential variations to circuit components with correlation constraints [8]. As electrical devices are typically small as compared to operation wavelengths, phase errors caused by variations are small; therefore, correlations are only applied to critical components that require matching, typically device pairs. For example, a differential amplifier requires matching resistor pairs, as its common-mode-rejection-ratio is sensitive to differential variations on the resistors [9]. However, in photonic designs, as devices are much longer than the operation wavelength, a small change in waveguide dimensions can lead to a dramatic phase error, especially for interferometric devices with long waveguide paths. Therefore, the MC approach requires to assign random variations differentially to all of the circuit components, and more importantly, all of the assigned variations are required to be correlated and layout-dependent. Not only does this approach face serious scaling issues (e.g., it requires N2 correlation parameters for N components), but it is also difficult to determine correlation parameters from a physical layout. Recently, photonics manufacturing variation analysis is driving more attention, and there have been models for dealing with the variability of SOI devices [10–12] and reports for the manufacturing correlations [13–16]. Nevertheless, photonic circuit simulation models that are capable of handling the layout-dependent correlated variations remain to be developed.
This paper addresses the two aforementioned challenges. In section 2, we propose a simple method to characterize manufacturing variations, based on the spectral response variations of a transverse electric (TE) mode racetrack resonator. Using the proposed method, we characterize the manufacturing variations of a 200-mm-wafer, which was fabricated using a 248-nm deep ultraviolet (DUV) lithography photonics process. In section 3, we develop models for our enhanced MC simulation methodology. Figure 1 sketches our simulation approach: a photonics layout is transferred to circuit simulator using netlist extraction [17]; based on characterization data from fabrication, wafer-scale correlated variations are simulated as virtual wafers on a discrete grid, and the variations are assigned to each circuit component according to the component’s layout coordinates obtained from the netlist extraction; each component updates its performance according to the obtained local physical dimensions, and therefore, circuit simulation can capture the layout-dependent correlated variations. The simulation flow and models will be detailed. In section 4, using the enhanced MC simulation, we present numerical experiments on several PIC designs to study their performance variability. Discussion and conclusion are given in sections 5 and 6, respectively.
Fig. 1 Proposed circuit simulation approach. First, the designer generates a layout and uses the netlist extraction tool to transfer the layout to a circuit simulator. Second, wafer-scale correlated variations are automatically simulated based on the characterized results from fabrication, and variations at each component position are assigned to the component. Then, the optical response of each component will automatically get updated according to obtained variations. Finally, circuit simulation is performed. Results include both common mode and differential mode variability.
2. Manufacturing variation characterization
2.1. Principle for characterization methodology
On SOI platforms, strip waveguides are widely used for on-chip optical routing due to their good confinement for optical propagation modes. However, the wavelength-dependent effective refractive index, neff (λ), of an optical mode in a strip waveguide is sensitive to the waveguide width and height. Using a mode solver, we calculated the neff (λ) for the fundamental transverse electric (TE) mode in strip waveguides with various widths and heights. The cross-section schematic for the simulated silicon waveguide is shown in Fig. 2(a), and simulation results of neff (λ) for various waveguide widths and heights are shown in Figs. 2(b) and 2(c), respectively.
Fig. 2 (a) Cross-section schematic of an SOI strip waveguide; (b) effective index of strip waveguides for various waveguide widths; (c) effective index of strip waveguides for various waveguide heights.
Racetrack resonators are interferometric devices, and manufacturing variations in the waveguide width, w, and height, h, both affect the neff (λ) of the waveguide, and consequently result in variations in the spectral response. In order to characterize manufacturing variations, an all-pass racetrack resonator was designed as a test device. The schematic layout for such a device is shown in Fig. 3(a), where silicon waveguides of the racetrack resonator have a nominal width of 500 nm and a nominal height of 220 nm. The radius, coupling length and coupling gaps of the device are 12 μm, 4.5 μm, and 200 nm, respectively. The through-port transmission spectrum of the device is shown in Fig. 3(b), and is given by [18]:
where t=0.9592 is the self-coupling coefficient at 1550 nm and is obtained using a finite-difference time-domain (FDTD) solver; L is the round trip length of the resonator; a=0.9811 is the round-trip amplitude attenuation factor contributed by both waveguide propagation loss in the resonator (assumed to be 10 dB/cm) and the loss in the coupler.Fig. 3 (a) Schematic layout of the racetrack resonator test device; (b) transmission spectrum for such a device without fabrication variations. FSR is free-spectral-range.
Manufacturing variability has a direct impact on the resonance wavelengths, λres, of racetrack resonators. As an example, Fig. 4(a) shows transmission spectra for racetrack resonators with various waveguide widths, from which resonance shifts can be directly observed. Manufacturing variability also affects the waveguide group index, ng, at each resonance wavelength λres. The ng can be extracted from the spectral response and is determined by [13, 19]:
where FSR is free-spectral-range and is defined as wavelength spans between two adjacent resonances, as illustrated in Fig. 3(b). By quantifying the variations of λres and ng versus the variations of waveguide width and height (i.e., , , and ), waveguide dimensions of a racetrack resonator can be extracted from its spectral response by: where Δw and Δh are deviations of waveguide width and height relative to their nominals, respectively; Δng and Δλres are variations of group index and resonance wavelength relative to their nominals at a resonance mode, respectively.Fig. 4 Simulation results for racetrack resonators with a nominal waveguide height of 220 nm and various waveguide widths. (a) Transmission spectra; (b) group indices at resonances; (c) resonance wavelength versus waveguide width for a selected resonance mode; (d) group index versus waveguide width for a selected resonance mode.
Using Eq. (2), we extract ng from the spectral responses of racetrack resonators with various waveguide widths, w, and a nominal waveguide height, h, of 220 nm. The extracted ng results are shown in Fig. 4(b). It is found that ng decreases with an increase of w whereas its corresponding λres increases. According to [13], the data points forming a downward diagonal line have the same resonance mode. Based on this relationship, it is possible to identify the data points that belong to the same mode and accordingly determine their resonance shifts even if the shifts are greater than the FSR. Here, we analyze a resonance mode having a nominal λres of 1544.25 nm and ng of 4.1899, as shown in 4(b). For this mode, we plot the λres versus w and ng versus w in Figs. 4(c) and 4(d), respectively, both of which can be approximated using linear fits. According to the fitting slopes, we obtain:
Similarly, we simulate racetrack resonators having a nominal waveguide width of 220 nm and various waveguide heights, and extract the ng at resonances from each simulated spectrum, as shown in Fig. 5(b). It is found that the data points for each resonance mode change in a upward diagonal direction, with an increase in h. For the resonance mode having a nominal λres of 1544.25 nm and a nominal ng of 4.1899, the λres versus h and ng versus h are plotted in Figs. 5(c) and 5(d), respectively. Both the λres versus h and ng versus h can be approximated using linear fits, and according to the fitting slopes we obtain:
Fig. 5 Simulation results for racetrack resonators with a nominal waveguide width of 500 nm and various waveguide heights. (a) Transmission spectra; (b) group indices at resonances; (c) resonance wavelength versus waveguide height for a selected resonance mode; (d) group index versus waveguide height for a selected resonance mode.
Based on Eq. (4) with data from Eqs. (5)–(8), the deviations of waveguide dimensions can be extracted from the variations of resonance wavelength and group index. While here the example test device for our characterization method is based on a racetrack resonator, ring resonators can also be used as test devices.
2.2. Characterization error
As Eqs. (5)–(8) are based on linear approximations, there are numerical errors in our proposed variation extraction method, and the errors can be quantified by using the following procedures. First, the transmission spectrum of a racetrack resonator is simulated with given waveguide width and height deviations, i.e., Δwgiven and Δhgiven, respectively. Second, the λres for each resonance mode is obtained from the simulated spectrum, and its corresponding ng is extracted by using Eq. (2). Third, the Δλres and Δng are calculated for the resonance mode having a nominal λres of 1544.25 nm and a nominal ng of 4.1899. Finally, base on Eq. (4) the deviations of waveguide width, Δwextracted, and waveguide height, Δhextracted, are extracted. The width extraction error, ErrorΔw, and height extraction error, ErrorΔh, are evaluated by:
For most of photonics 193-nm and 248-nm deep ultraviolet (DUV) lithography processes, the etching linewidth is typically controlled to an accuracy of within ±20 nm, and the wafer thickness uniformity is guaranteed to an accuracy of within ±10 nm. In our error tests, extraction errors are investigated in the same ranges. Figures 6(a) and 6(b) show the extraction ratios for Δw and Δh, respectively. It is found that the ErrorΔw is less than 0.85 nm for a maximum Δwgiven of 20 nm and the ErrorΔh is less than 0.55 nm for a maximum Δhgiven of 10 nm.
Fig. 6 Error test results in a Δw deviation range of ±20 nm and a Δh deviation range of ±10 nm, for the proposed variability characterization method. (a) Width extraction error; (b) height extraction error.
2.3. Experiment results
Using the proposed methodology, we characterized the manufacturing variability of a 200-mm-wafer fabricated by a 248-nm DUV lithography process through IME’s silicon photonics foundry. A total number of 2074 identical racetrack resonators were fabricated on 34 lithography dies on the wafer. The size for each die is 7.6 mm × 6 mm, and each die has 61 identical devices. Figure 7(a) shows the schematic layout for a fabricated racetrack resonator, the design parameters of which are the same as those of the design as shown in Fig. 3(a). Figures 7(b) and 7(c) show the layout distributions for the racetrack resonators on each die and the wafer map for the fabricated 200-mm-wafer, respectively.
Fig. 7 (a) Schematic layout for the fabricated racetrack resonator; (b) distribution of racetrack resonators on each wafer die; (c) wafer map for the fabricated multi-project-wafer.
All of the fabricated racetrack resonators were measured using an automated photonics testing setup, and the temperature of the setup was stabilized with a control stability of 1 mK. Figure 8(a) shows the spectra of devices on die #20, which is located close to the centre of the wafer. Using Eq. (2), we extracted the ng of the devices on die #20, and results are shown in Fig. 8(b). Each cluster of data points in Fig. 8(b) represents a resonance mode. The data points for the mode being closest to 1544.25 nm are selected for Δw and Δh extractions. In the extractions, first, the Δλres and Δng for each selected data point are calculated relative to a nominal λres of 1544.25 nm and a nominal ng of 4.1899 for the selected mode. Then, Δw and Δh are obtained using Eq. (4). Figures 8(c) and 8(d) respectively show the extracted Δw and Δh versus coordinates on the die #20. It is found that Δw varies from 3.15 nm to 9.14 nm; while Δh varies from −0.44 nm to −3.03 nm.
Fig. 8 Characterization results for the manufacturing variability a 200-mm-wafer. (a) Measured spectra for the 61 devices on die #20; (b) extracted ng for the 61 devices on die #20; (c) extracted Δw versus position on die #20; (d) extracted Δh versus position on die #20; (e) extracted ng versus λres for the 2074 devices on the wafer; (f) histogram for extracted Δw across the wafer; (g) histogram for extracted Δh across the wafer.
Figure 8(e) shows the extracted ng versus λres for all of the 2074 devices on the wafer, and each cluster of data points corresponds to a resonance mode. As the deviations of λres and ng in the wafer-scale are large, it is difficult to determine the modes for the data points being close to mode boundaries. However, by shrinking the size of data points selection region by 20% as illustrated by red dash lines in Fig. 8(e), it is more likely that the data points within the 80% region belong to the same mode. Accordingly, we select the data points within the 80% region for Δw and Δh extractions, in which the nominal λres and ng for the mode of the selected data points are 1544.25 nm and 4.1899, respectively. The statistical results for extracted Δw and Δh are shown in Figs. 8(f) and 8(g), respectively. We fit the histograms using normal distributions to estimate the variations across the 200-mm-wafer, and fitting results are listed in Table 1. According to the fitting results, the standard deviations for waveguide widths and heights are 3.855 nm and 1.316 nm, respectively.
Table 1. Statistical results for the manufacturing variations of a 200-mm-wafer fabricated through a 248-nm DUV lithography process
Having more test devices on the wafer can increase the amount of data samples and consequently improve the accuracy of histogram fitting results. In addition, using test devices with a larger FSR [20] can avoid the mode overlapping issue, and as a result, make it easy to determine the resonance mode of each extracted ng data points.
2.4. Statistical results for manufacturing variations in literature
In recent years, there has been studies investigating the manufacturing variations [3,4,6,21–25] and correlations [13, 14] of various photonics fabrication processes. In these reports, many techniques including SEM, ellipsometry, AFM, and indirect measurements were used for characterizations. Here, we summarize the characterized wafer-to-wafer and within-wafer variations in Tables 2 and 3, respectively. Wafer-to-wafer variations evaluate the drifts for the mean of waveguide linewidth and height. As demonstrated in [4], the standard deviation for the mean of waveguide height from wafer to wafer is 1.26 nm. For within-wafer variations, the standard deviation for waveguide linewidth, σΔw, ranges from 0.78 nm to 2.65 nm; while that for waveguide height, σΔh, ranges from 0.83 nm to 4.16 nm. According to [13, 14], within-wafer variations are correlated over short distance scales, and the characterized correlation length, l, is around 4 mm to 5 mm.
Table 2. Literature results for wafer-to-wafer fabrication variations
Table 3. Literature results for within-wafer fabrication variations
3. layout-dependent, enhanced Monte Carlo simulation methodology
As discussed in the introduction, photonics circuit simulations need to take layout-dependent correlated variations into account in order to predict the reliability of circuit performance. This section discusses simulation flows and models for such circuit simulations.
3.1. Simulation procedures
The principle idea of our enhanced MC simulation is simulating spatially correlated manufacturing variations in a wafer scale and mapping the variations to each circuit component according to its layout coordinate, so that circuit simulations can capture the correlated variations between devices. Figure 9(a) shows the flow chart of the proposed methodology. The procedures are as follows:
- Layout to schematic: First, the layout of a silicon photonic circuit under test is created (in the example presented here, it is designed using KLayout [27], an open-source layout tool). The netlist of the layout, which provides coordinates for primitive elements in the circuit, connections between the elements, and path and dimension information of waveguides, is extracted from the layout using a developed open-source tool [17, 28], and is then imported into a circuit simulator (in the example presented here, we use Lumerical INTERCONNECT [29]). As an example, Figs. 10(a) and 10(b) illustrate the physical layout and the simulation schematic of a simple photonic circuit consisting of two grating couplers connected by a strip waveguide.
- Manufacturing variation simulation: Second, in the circuit simulator, we create a virtual wafer for waveguide width variations, Δw, and a virtual wafer for waveguide height variations, Δh, as illustrated in Figs. 9(b) and 9(c), respectively. Each virtual wafer is characterized by a standard deviation amplitude and a correlated length, both of which are based on experimental results [3, 4, 6, 13, 14, 21–25]. The virtual wafer model will be described in section 3.2.
- Variation mapping: The simulated virtual wafers, as illustrated in Figs. 9(b) and 9(c), can be divided into numerous lithography dies with identical size. The virtual maps on one of the wafer dies is selected and mapped to the photonic circuit under test. As an example, Fig. 11 illustrates the die selection and variation mapping for a ring resonator circuit. Each circuit component obtains its local Δw and Δh according to its layout coordinates.
- Performance interpolation: The performance of each circuit component is updated according to the obtained Δw and Δh, which is performed using interpolation in the circuit simulator. The parameterized component models will be described in section 3.3.
- Circuit simulation: The photonic circuit, with updated components’ performance, is simulated, and simulation results can be visualized.
Fig. 9 (a) Flow chart for the enhanced Monte Carlo methodology; (b) simulated virtual wafer for waveguide width deviations; (c) simulated virtual wafer for waveguide height deviations.
Fig. 10 An example circuit consisting of two grating couplers connected by a waveguide. (a) Physical layout; (b) simulation schematic in the circuit simulator.
The proposed methodology provides a within-wafer analysis and a wafer-to-wafer analysis. The within-wafer analysis iterates the steps of “variation mapping”, “performance interpolation”, and “circuit simulation”, to study the circuit performance variations across various lithography dies, as shown in Fig. 9(a). Figure 11 illustrates the die selection and variation mapping for the within-wafer analysis. The wafer-to-wafer analysis includes the within-wafer analysis, and each of its iteration generates a new Δw virtual wafer and a new Δh virtual wafer for the within-wafer analysis, as illustrated in Fig. 12. The iteration numbers for the within-wafer analysis and wafer-to-wafer analysis, i.e., N and M respectively as shown in Fig. 9(a), can be defined by users.
3.2. Virtual wafer model
We model the within-wafer manufacturing variations using a correlated surface roughness function [30]. Firstly, we generate an uncorrelated random distribution map, z(x, y), based on discrete mesh of points in x–y plane, and the value at each discrete point is a random number following a normal distribution with a mean of 0 and a standard deviation of σ. Then, the random height function z(x, y) is convolved with a Gaussian filter, g(x, y), which is given by:
in which, l is correlation length for variations. The convolved wafer with correlated variation is given by: where ℱ and ℱ− are denotations for the fast Fourier transform and the inverse fast Fourier transform, respectively. As an example, Figs. 13(a)–13(c) show a random distribution map z(x, y), a Gaussian filter map g(x, y), and a correlated variation wafer map m(x, y), respectively.Fig. 13 (a) A 100 mm × 100 mm random distribution map z(x, y) generated with σ = 2; (b) a Gaussian filter map g(x, y) generated with l = 4 mm; (c) the correlated variation wafer map m(x, y).
We include wafer-to-wafer variations in our wafer model by adjusting the mean of the correlated variation wafer map, which is given by:
where is c a random number following a normally distribution with a mean of 0 and a standard deviation σ∗. In total, there are three input parameters to our virtual wafer model, which are σ, l, and σ∗.To increase computation efficiency, virtual wafers are simulated using a coarse simulation mesh of 500 μm × 500 μm in the “manufacturing variation simulation” step. In the “variation mapping” step, the simulated width and height variations on the selected die are interpolated to have a high resolution mesh of 5 μm × 5 μm. The interpolated high resolution variation maps are then mapped to photonic circuits. As an example, Fig. 14 illustrates the variations maps before and after interpolation.
Fig. 14 (a) A 100 mm × 100 mm correlated variation wafer map m(x, y), which is simulated using a coarse simulation mesh; (b) and (c) are the variation maps for a 10 mm × 10 mm die located at the top right corner of the wafer before and after interpolation, respectively.
3.3. Parameterized component models
In order to consider fabrication variations in circuit simulations, we require parameterized component models that are continuous with width and height variations. In addition, wavelength-dependency is necessary in order to model a circuit over a wavelength range of interest. This section describes how waveguides and components are modeled.
3.3.1. Waveguide compact model
Waveguides are fundamental building blocks in PICs. In a circuit simulator, such as Lumerical INTERCONNECT [29], waveguide properties can be described using waveguide loss and wavelength-dependent effective refractive index, neff (λ). In our waveguide model, the waveguide loss is based on empirical results (e.g. 2–6 dB/cm for 500 nm ×220 nm waveguides), and it can be defined by the user. The neff (λ) can be described using a third-order Taylor expansion, which is given by [19]:
where λ0 = 1550 nm is the central wavelength of operation bandwidth; and are determined by the group index ng and group velocity dispersion D [19], respectively:The neff (λ0), ng (λ0) and D(λ0) can be obtained from eigenmode simulations.
To parameterize the neff (λ) as a function of the waveguide width w and waveguide height h, we perform eigenmode simulations for waveguides with various geometries (e.g. sweeping w from 300 nm to 1000 nm with a 20 nm step, and sweeping h from 200 nm to 240 nm with a 2 nm step), and we generate a look-up table for the simulated neff (λ0), ng (λ0) and D(λ0) of each waveguide geometry. According to the look-up table, we use a multi-dimensional spline interpolation method to interpolate the neff (λ0), ng (λ0) and D(λ0) for the desired w and h, and accordingly obtain the neff (λ) using Eq. (13). The look-up table and the interpolation are implemented within our waveguide compact model in Lumerical INTERCONNECT.
In our waveguide model, the impacts of the spatially dependent Δw and Δh on a waveguide are averaged along its layout path. The waveguide path information is extracted from the layout, as discussed in section 3.1, and is input into our waveguide model. As illustrated in Fig. 15, with the path information, a waveguide is mathematically meshed into pieces by the simulated virtual wafers, and each waveguide piece automatically obtains the Δw and Δh at its mesh point. The average waveguide width and height along the path, i.e., wavg and havg, are given by:
where N is total number of waveguide pieces; Δwn and Δhn are width and height variations for the waveguide section n, respectively. The averaging is implemented as a script within our waveguide model. Based on wavg and havg, the neff (λ0), ng (λ0) and D(λ0) are determined as described above.3.3.2. S-parameter component models
All of the primitive components, such as y-junction splitter [31], grating couplers [32], and directional couplers [33], can be described using optical S-parameters [19, 29], which can be calculated by using commercially available Finite-Difference Time-Domain (FDTD) solvers, such as Lumerical FDTD Solutions [29]. Optical S-parameters include the amplitude and phase responses of a primitive component (indicating both transmission and reflection parameters). Since these simulations do not provide propagation loss result from waveguide sidewall roughness, we need to add this separately to the compact model. To parameterize the primitive component model as a function of Δw and Δh, first, we obtain S-parameters for the nominal design and the four process corners of manufacturing variations. We consider a maximum Δw of ±20 nm and a maximum Δh of ±10 nm for a photonics process. Accordingly, the four corners for (Δw, Δh) are (20 nm, 10 nm), (−20 nm, 10 nm), (−20 nm, −10 nm) and (20 nm, −10 nm), and the nominal design is (0 nm, 0 nm), as illustrated in Fig. 16(a). Then, for a desired (Δw, Δh) input, we use a multi-dimensional spline interpolation method to interpolate the amplitude and phase of the S-parameters. After interpolation, the S-parameters’ passivity and reciprocity are tested and enforced [19], to avoid violating energy conservation laws. We implement the S-parameter data and interpolation as a script within the component model in Lumerical INTERCONNECT. As an example, Fig. 16(b) shows the S31 amplitudes for the nominal design, four process corners and one interpolated case for a TE mode 2×2 directional coupler. The example directional coupler has a coupling length of 17.5 μm, a coupling gap of 200 nm, waveguide widths of 500 nm, and waveguide heights of 220 nm.
Fig. 16 S-parameter component model. (a) Process corners for two process parameters: waveguide width variation, Δw, and height variation, Δh; (b) S-parameters for the nominal design, four process corners, and one interpolated geometry of a 2×2 splitter.
3.3.3. Sub-circuit components
Sub-circuit components can be described using both the S-parameter model and the waveguide model. One example of a sub-circuit component is a ring resonator. An all-pass ring resonator can be decomposed into a directional coupler described by the S-parameter model and a bent waveguide described by the waveguide model, as shown in Fig. 17(a); An add-drop ring resonator can be decomposed into two directional couplers both described by S-parameter models, as shown in Fig. 17(b). The performance of the ring resonator, such as resonance wavelength, FSR and extinction ratio, will vary with the thickness and height variations introduced into the directional coupler model and the waveguide model. Numerical experiments for ring resonators will be presented in section 4.
Fig. 17 Layout decomposition examples of ring resonators. (a) All-pass ring resonator; (b) add-drop ring resonator.
4. Numerical experiments
In this section, we study the performance variations of several PIC designs within a single 200 mm × 200 mm wafer by using the proposed enhanced MC simulation methodology that is described in section 3. The input parameters to the virtual wafer model of MC simulations are listed in Table 4. The wafer size and die size in our numerical experiments are 200 mm × 200 mm and 5 mm × 5 mm, respectively. In total, the wafer includes 1600 dies. Standard deviations for the waveguide width and height, i.e., σΔw and σΔh, are set to 3.855 nm and 1.316 nm, respectively, which are based on the measurement results of the 200-mm-wafer as presented in section 2.3. The correlation lengths for the width and height variations are both 4.5 mm, according to [13, 14]. Note that simulation results of these examples are dependent on the input parameters to the virtual wafer model, and these parameters can be user-defined according to the variability of a specific fabrication process.
Table 4. Input parameters for the virtual wafer model
4.1. Ring resonators based single-channel filter
Our first numerical experiment investigates a ring resonator based single-channel filter on the SOI platform with a 220 nm thick silicon layer. Figure 18(a) shows the layout of the filter design, which is generated using Klayout [27]. The filter circuit includes an add-drop ring resonator, three grating couplers [32] as inputs and outputs, and silicon waveguides connecting the ports of the resonator and those of the grating couplers. One of the resonator ports is terminated using a tapered waveguide terminator. The ring radius, waveguide width, and coupling gaps of the ring resonator are designed to be 5 μm, 500 nm, and 100 nm, respectively. All of the connecting waveguides are 500 nm wide. Using an open-source netlist extraction tool [28] that was detailed in our previous work [17], we extract the netlist of the filter circuit from its physical layout, and load it into Lumerical INTERCONNECT for circuit simulations. Figure 18(b) shows the loaded simulation schematic of the filter circuit. All of the circuit components have been parameterized using the methods as described in section 3.3.
Fig. 18 A ring resonator based single-channel photonic filter circuit. (a) Physical layout; (b) circuit simulation schematic.
Figure 19(a) shows the simulated transmission spectra at the through-port and the drop-port of the filter, without having fabrication variation. Ideally, the filter has a resonance wavelength λres at 1542.07 nm, with an extinction ratio (ER) of 17.8 dB. We define ER as the absolute difference, in a logarithm scale, between the through port transmission and drop port transmission at a resonance.
Fig. 19 Simulation results for a ring resonator based single-channel filter. (a) ideal transmission spectrum; (b) transmission spectra on 3 representative dies with fabrication variations; (c) extracted group index versus resonance wavelength for the filter in 1600 simulations; (d) and (e) are histograms of resonance wavelengths and extinction ratios for the filter in 1600 simulations, respectively.
Manufacturing variations affect the filter performance in numerous ways. Grating couplers are sensitive to the width and height variations on the grating teeth, both of which vary their insertion losses and central wavelengths. For the ring resonator, the waveguide dimension variations change its coupling strengths, and therefore, result to variations on the ERs. The waveguide dimension variations also vary λres. We simulate the filter performance on all of the 1600 dies in the wafer, and obtain through-port and drop-port transmission spectra from each circuit simulation. Figure 19(b) shows transmission spectra of the filter on 3 representative dies, where resonance shift, ER variation, and insertion loss variation are observed. To statistically analyze the filter performance, we extracted group indices ng from all of the 1600 simulated spectra by using Eq. (2), and the results are shown in Fig. 19(c). Each cluster of data points in Fig. 19(c) represents a resonance mode, so we can analyze the performance variations of the filter based on the data points of the same resonance mode. Figures 19(d) and 19(e) show histograms for the λres and ERs of a resonance mode, respectively. Normal distribution fits are used to analyze the histograms, as shown by the red lines in Figs. 19(d) and 19(e), and fitting results are summarized in Table 5. According to Table 5, the predicted mean and standard deviation for λres are 1541.9 nm and 2.799 nm, respectively; the predicted mean and standard deviation for ERs are 17.65 dB and 0.555 dB, respectively.
Table 5. Monte Carlo analysis results for the ring resonator based single-channel filter
Based on the results of MC analysis, designers can develop appropriate strategies to compensate the fabrication variations. For example, a thermal tuner can be used to adjust the λres of the filter, and the maximum power consumption of the thermal tuner can be estimated according to the shifting range of λres.
4.2. Ring resonator based dual-channel filter
As demonstrated in many reports [13, 14], fabrication variations are correlated over short distance scales, and as a result, adjacent devices in layouts tend to have similar performance variations. One example to show the correlated variations is a dual-channel filter based on two ring resonators. Figure 20(a) shows the layout for such a device, which consists of two add-drop ring resonators, grating couplers [32] as inputs and outputs, waveguide terminators, and connecting waveguides. Both rings have waveguide widths of 500 nm and coupling gaps of 100 nm. The two rings are designed to have a different radius, which is 5 μm for the ring 1 and 5.05 μm for the ring 2, so that they have different resonance channels. The center-to-center distance, d, between the two rings in the layout is 30 μm, as shown in Fig. 20(a).
Fig. 20 A photonic ring resonator based dual-channel filter circuit. (a) Physical layout; (b) circuit simulation schematic.
First, we perform a circuit simulation without introducing fabrication variations. The simulated spectra at both the through port and the drop ports of the filter are shown in Fig. 21(a). Ideally, the designed λres of Channel 1 and Channel 2 are 1542.07 nm and 1551.06 nm, respectively. We define channel spacing, Δλ, as the difference between the two resonance wavelengths, which is 9.01 nm. Next, we simulate the filter performance on all of the 1600 dies of the wafer, and obtain output transmission spectra from each simulation. Statistical results for the λres of the two channels and the Δλ between the two channels are shown in Figs. 21(b) and 21(c), respectively. According to Fig. 21(b), Channel 1 shifts in between 1533 nm to 1552 nm; while Channel 2 shifts in between 1541 nm to 1561 nm. However, Δλ only varies in a small range from 8.8 nm to 9.15 nm, as shown in Fig. 21(c), which indicates that the shifts of the two channels are strongly correlated. Table 6 summarizes the parameters for the histogram fits of Figs. 21(b) and 21(c).
Fig. 21 Simulation results for the ring resonator based dual-channel filter. (a) Output transmission spectra without fabrication variation; (b) and (c), respectively, are histograms for resonance wavelengths λres and channel spacings Δλ of the filter in 1600 simulations.
Table 6. Monte Carlo analysis results for a dual-channel filter design with d=30 μm.
Another interesting numerical experiment is to investigate the dependence of correlated variations on the spatial distances between devices. Here, we study the variations of Δλ for filter designs with various center-to-center distances, which are d=30 μm, d=150 μm, and d=300 μm. For each filter design, we simulate the filter performance on all of the 1600 dies on the wafer, and analyze the simulated results. Figure 22 shows and compares the variations of Δλ for the three designs. We fit the histograms using normal distributions, and the fitting results are listed in Table 7. According to the results, a larger center-to-center distance results in a larger variation on the channel spacing because fabrication variations of two devices are less correlated when the devices are further apart. This numerical experiment indicates the significance of a compact layout for PIC designs, when performance matching between components is necessary.
Fig. 22 Histograms for the channel spacings for filter designs with d=30 μm, d=150 μm, and d=300 μm.
Table 7. Monte Carlo analysis results for various dual-channel filters designs
5. Discussion
In electronic digital circuits, metrics for circuit yield can be based on the yield of each digital component. Typically, a digital circuit will work if all of its components operate properly. However, most photonic designs are analog circuits; as a result, photonics yield analysis is based on circuit performance rather than single device performance, which is similar to that of analog electronics. For example, when considering the photonic two-channel filter that is presented in section 4.2, designers need to define circuit performance metrics for yield (e.g., channel spacing of the two channels being 9 nm +/− 0.1 nm), rather than simply evaluating the performance of a single component (e.g., wavelength variations of a single channel). Similarly, for a transceiver (transmitter and receiver), the performance metric may be bit error rate (BER); in this case, a threshold BER can be used to determine yield.
In photonics fabrications, local pattern density can affect both lithography and etching processes, and therefore, it results in devices with different pattern density having different variations even if they are adjacent in the layout, e.g., a fabricated SOI arrayed waveguide grating (AWG) may have different variability from a ring resonator. This effect is not considered in our proposed model but can be included in our future work.
One possible solution to deal with local pattern density effects is to treat pattern density on a circuit basis. A location-dependent pattern density map of a photonic circuit can be obtained by scanning a sampling window across the layout and calculating the layer density at each sampling spot (2D low-pass filter, with an appropriate parameters determined by characterizing the fabrication process). As an example, Fig. 23(a) shows the layout of a Mach-Zehnder interferometer (MZI) circuit, where one of its arms is surrounded by silicon tiling structures while the other is not. Figure 23(b) shows the calculated location-dependent pattern density map of the MZI circuit. It can be seen that the two arms have different pattern density along their paths, and therefore may have different resulting fabricated feature sizes. The calculated location-dependent pattern density can be transformed into waveguide geometric variations based on known parameters for the pattern density effects (which needs to be provided by the foundries; we note that silicon photonics foundries currently do not provide this information in their Process Design Kits, hence there is substantial work required to determine the fabrication process parameters relating to local pattern density). The circuit simulation would proceed as proposed here, where the performance of each circuit component is interpolated based on local variations from both the local pattern density effects and the Monte Carlo virtual wafers.
Fig. 23 An example Mach-Zehnder interferometer circuit. (a) Circuit layout; (b) local pattern density for the Si layer of the layout, which is calculated using a 6 μm × 6 μm sampling window with a scanning step of 0.3 μm.
6. Conclusion
This paper addressed two challenges of PIC designs: characterizing manufacturing variations and predicting circuit performance under the impacts of correlated manufacturing variations. We have demonstrated a simple and accurate characterization method for the manufacturing variations of PICs, which extracts the waveguide width and height from the spectral response of a single racetrack resonator. Using this method, we have measured the spectral responses of 2074 identical racetrack resonators on a 200-mm-wafer that was fabricated through a 248-nm DUV lithography photonics process, and we extracted the waveguide width and height at each device location on the wafer. According to the statistical results, the standard deviations for waveguide widths and heights for such a fabrication process are 3.855 nm and 1.316 nm, respectively. In addition, we have developed simulation flows and models to predict the performance of PICs under the impact of correlated manufacturing variations. The principle idea of our approach is modelling spatially-correlated variations on a discrete grid and assigning the variations to each circuit component so that circuit simulations can capture layout-dependent correlated variations. The simulation flows and theoretical models for the enhanced MC simulation methodology have been detailed. Using the enhanced MC simulation, we have presented numerical experiments on several PIC designs to study their performance variability.
To facilitate the proposed MC simulation, we have implemented a design and analysis toolbox that is integrated into an open-source layout tool, KLayout [27]. This toolbox is open-source and is available for download [28]. This toolbox offers a layout-first design methodology. It includes netlist extraction routines to generate a schematic from the layout and directly sends the netlist information to a circuit simulator, allowing users to perform MC simulations from the layout with user-defined manufacturing variations.
Acknowledgments
We acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC), including the Silicon Electronic-Photonic Integrated Circuits (SiEPIC) Program, and acknowledge CMC Microsystems for fabrication access and design tools. We appreciate discussions with Mr. Chris Cone and Dr. John Ferguson from Mentor Graphics, Dr. Stefan Preble from Rochester Institute of Technology, and Dr. Dan Deptuck from CMC Microsystems.
References and links
1. W. Bogaerts, M. Fiers, and P. Dumon, “Design challenges in silicon photonics,” IEEE J. Sel. Top. Quantum Electron. 20, 1–8 (2014). [CrossRef]
2. L. Chrostowski, Z. Lu, J. Flückiger, J. Pond, J. Klein, X. Wang, S. Li, W. Tai, E. Y. Hsu, C. Kim, J. Ferguson, and C. Cone, “Schematic driven silicon photonics design,” Proc. SPIE 9751, 975103 (2016). [CrossRef]
3. S. K. Selvaraja, W. Bogaerts, P. Dumon, D. V. Thourhout, and R. Baets, “Subnanometer linewidth uniformity in silicon nanophotonic waveguide devices using CMOS fabrication technology,” IEEE J. Sel. Top. Quantum Electron. 16, 316–324 (2010). [CrossRef]
4. S. K. Selvaraja, “Wafer scale fabrication technology for silicon photonic integrated circuit,” Ph.D. thesis, Ghent University (2011).
5. W. A. Zortman, D. C. Trotter, and M. R. Watts, “Silicon photonics manufacturing,” Opt. Express 18, 23598–23607 (2010). [CrossRef] [PubMed]
6. N. Ayotte, A. D. Simard, and S. LaRochelle, “Long integrated Bragg gratings for SOI wafer metrology,” IEEE Photonics Technol. Lett. 27, 755–758 (2015). [CrossRef]
7. X. Wang, W. Shi, H. Yun, S. Grist, N. A. F. Jaeger, and L. Chrostowski, “Narrow-band waveguide Bragg gratings on SOI wafers with CMOS-compatible fabrication process,” Opt. Express 20, 15547–15558 (2012). [CrossRef] [PubMed]
8. L. Lavagno, L. Scheffer, and G. Martin, EDA for IC Implementation, Circuit Design, and Process Technology (CRC, 2006).
9. A. S. Sedra and K. C. Smith, Microelectronic Circuits (Oxford University, 1998).
10. 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). [CrossRef] [PubMed]
11. D. Melati, A. Melloni, and F. Morichetti, “Real photonic waveguides: guiding light through imperfections,” Adv. Opt. Photonics 6, 156–224 (2014). [CrossRef]
12. 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). [CrossRef]
13. L. Chrostowski, X. Wang, J. Flueckiger, Y. Wu, Y. Wang, and S. T. Fard, “Impact of fabrication non-uniformity on chip-scale silicon photonic integrated circuits,” in Optical Fiber Communication Conference (2014), paper. Th2A.37. [CrossRef]
14. Y. Yang, Y. Ma, H. Guan, Y. Liu, S. Danziger, S. Ocheltree, K. Bergman, T. Baehr-Jones, and M. Hochberg, “Phase coherence length in silicon photonic platform,” Opt. Express 23, 16890–16902 (2015). [CrossRef] [PubMed]
15. A. Agarwal, D. Blaauw, and V. Zolotov, “Statistical timing analysis for intra-die process variations with spatial correlations,” in International Conference on Computer Aided Design (2003), pp. 900–907.
16. R. Wu, C. H. Chen, T. C. Huang, R. Beausoleil, and K. T. Cheng, “Spatial pattern analysis of process variations in silicon microring modulators,” in IEEE Optical Interconnects Conference (2016), pp. 116–117.
17. L. Chrostowski, Z. Lu, J. Flueckiger, X. Wang, J. Klein, A. Liu, J. Jhoja, and J. Pond, “Design and simulation of silicon photonic schematics and layouts,” Proc. SPIE 9891, 989114 (2016). [CrossRef]
18. W. Bogaerts, P. De Heyn, T. Van Vaerenbergh, K. De Vos, S. Kumar Selvaraja, T. Claes, P. Dumon, P. Bienstman, D. Van Thourhout, and R. Baets, “Silicon microring resonators,” Laser Photonics Rev. 6, 47–73 (2012). [CrossRef]
19. L. Chrostowski and M. Hochberg, Silicon Photonics Design (Cambridge University, 2015). [CrossRef]
20. N. Eid, H. Jayatilleka, M. Caverley, S. Shekhar, L. Chrostowski, and N. A. F. Jaeger, “Wide FSR silicon-on-insulator microring resonator with bent couplers,” IEEE 12th International Conference on Group IV Photonics (GFP, 2015), paper 96–97.
21. S. K. Selvaraja, E. Rosseel, L. Fernandez, M. Tabat, W. Bogaerts, J. Hautala, and P. Absil, “SOI thickness uniformity improvement using wafer-scale corrective etching for silicon nano-photonic device,” in Proceedings of the 2011 Annual Symposium of the IEEE Photonics Benelux Chapter (2011), pp. 289–292.
22. S. K. Selvaraja, E. Rosseel, L. Fernandez, M. Tabat, W. Bogaerts, J. Hautala, and P. Absil, “SOI thickness uniformity improvement using corrective etching for silicon nano-photonic device,” in 8th IEEE International Conference on Group IV Photonics (2011), pp. 71–73.
23. X. Chen, M. Mohamed, Z. Li, L. Shang, and A. R. Mickelson, “Process variation in silicon photonic devices,” Appl. Opt. 52, 7638–7647 (2013). [CrossRef] [PubMed]
24. S. K. Selvaraja, G. Winroth, S. Locorotondo, G. Murdoch, A. Milenin, C. Delvaux, P. Ong, S. Pathak, W. Xie, G. Sterckx, G. Lepage, D. Van Thourhout, W. Bogaerts, J. Van Campenhout, and P. Absil, “193nm immersion lithography for high-performance silicon photonic circuits,” Proc. SPIE 9052, 90520F (2014).
25. R. G. Beausoleil, A. Faraon, D. Fattal, M. Fiorentino, Z. Peng, and C. Santori, “Devices and architectures for large-scale integrated silicon photonics circuits,” Proc. SPIE 7942, 794204 (2011). [CrossRef]
26. D. X. Xu, J. H. Schmid, G. T. Reed, G. Z. Mashanovich, D. J. Thomson, M. Nedeljkovic, X. Chen, D. V. Thourhout, S. Keyvaninia, and S. K. Selvaraja, “Silicon photonic integration platform–have we found the sweet spot,” IEEE J. Sel. Top. Quantum Electron. 20, 189–205 (2014). [CrossRef]
28. https://github.com/lukasc-ubc/SiEPIC_EBeam_PDK
30. R. L. Wagner, J. Song, and W. C. Chew, “Monte Carlo simulation of electromagnetic scattering from two-dimensional random rough surfaces,” IEEE Trans. Antennas Prop. 45, 235–245 (1997). [CrossRef]
31. Y. Zhang, S. Yang, A. E.-J. Lim, G.-Q. Lo, C. Galland, T. Baehr-Jones, and M. Hochberg, “A compact and low loss y-junction for submicron silicon waveguide,” Opt. Express 21, 1310–1316 (2013). [CrossRef] [PubMed]
32. Y. Wang, X. Wang, J. Flueckiger, H. Yun, W. Shi, R. Bojko, N. A. F. Jaeger, and L. Chrostowski, “Focusing sub-wavelength grating couplers with low back reflections for rapid prototyping of silicon photonic circuits,” Opt. Express 22, 20652–20662 (2014). [CrossRef]
33. Z. Lu, H. Yun, Y. Wang, Z. Chen, F. Zhang, N. A. F. Jaeger, and L. Chrostowski, “Broadband silicon photonic directional coupler using asymmetric-waveguide based phase control,” Opt. Express 23, 3795–3808 (2015). [CrossRef] [PubMed]
