Accurate extraction of fabricated geometry using optical measurement

Yufei Xing; Jiaxing Dong; Sarvagya Dwivedi; Umar Khan; Wim Bogaerts

doi:10.1364/PRJ.6.001008

Abstract

We experimentally demonstrate extraction of silicon waveguide geometry with subnanometer accuracy using optical measurements. Effective and group indices of silicon-on-insulator (SOI) waveguides are extracted from the optical measurements. An accurate model linking the geometry of an SOI waveguide to its effective and group indices is used to extract the linewidths and thicknesses within respective errors of 0.37 and 0.26 nm on a die fabricated by IMEC multiproject wafer services. A detailed analysis of the setting of the bounds for the effective and group indices is presented to get the right extraction with improved accuracy.

1. INTRODUCTION

The submicron silicon-on-insulator (SOI) platform for silicon photonics offers tight confinement of light and compact integration of photonic devices. However, the high material index contrast also makes devices very sensitive to the geometry variation ^[1]. The variation introduced in fabrication often significantly deteriorates the device performance. Especially for spectral filters, geometry variation causes a shift in the spectrum and needs good compensation ^[2,3]. Therefore, an accurate evaluation of the fabricated geometry helps to estimate how to compensate the performance error and make a sensible design.

Extracting the fabricated linewidth and layer thickness is essential in getting the input data for analysis, such as performance evaluation ^[4,5] variability analysis ^[6,7], and revising compact behavior models ^[8]. However, metrology measurement of a fabricated photonic chip using a scanning electron microscope (SEM) is both expensive and destructive. Conventionally, the semiconductor fab only takes a few destructive cross-section pictures at given wafer locations on different dies, and this is only during process development, not in production. The geometry accuracy of such SEM-based measurements is usually good for a qualitative confirmation but not accurate enough for exact modeling. Ultraprecise methods such as the atomic force microscope (AFM) are extremely time-consuming ^[9]. In production, nondestructive methods such as top-down SEM, ellipsometry, and scatterometry are used to measure geometries such as linewidth and layer thickness. The problem with such measurements is that they are collected in dedicated measurement sites, which are often not representative of the actual waveguide devices. Current methods are not capable of extracting accurate wafer maps of the actual device geometry on a nanometer scale and its variability.

An alternative approach is to use optical transmission measurements on the actual devices to extract geometry parameters. Investigating the spectral response of devices such as microdisks or long Bragg gratings offers a more efficient way of characterizing manufacturing variations. However, most demonstrated methods either use dual-polarization measurements or request complex spectral reflection measurements from both ends of the device ^[4,10]. Optical properties such as the effective index and group index can be extracted from interfering structures such as Mach–Zehnder interferometers (MZIs) ^[11] and ring resonators ^[5,12]. Recent research shows the possibility to correlate waveguide geometry with these behavioral parameters such as resonance wavelength and group index by mapping them with a linear model ^[5]. Lu et al. measured spectral responses of ring resonators over the wafer and from this derived a geometry wafer map, demonstrating the potential of the optical method for wafer-scale geometry extraction. In this study, they used a ring resonator that consists of both straight and bent waveguides. Since straight and bent waveguides have different effective and group indices, the geometrical cross section of a straight waveguide cannot be extracted accurately from a ring, without making assumptions on the straight and bent geometry that cannot be verified.

Sign up for Photonics Research TOC. Get the latest issue of Photonics Research delivered right to you！Sign up now

Even though geometry extraction through transmission measurements shows real practical potential, it is still challenging to characterize the manufacturing quality accurately. First of all, the geometry model that links geometry with the behavior parameters should on one hand be very accurate, but on the other hand have as few parameters as possible. This is hard to achieve with combinations of different waveguide types (straight and bent), each with their own optical properties. Second, there is always noise in the spectral measurements, error in measurement alignment, and variability in grating couplers (GCs). It is essential to develop a method that is tolerant of the above-mentioned factors. Third, geometry extraction requires a method to simultaneously extract both effective and group indices from the same device. Most importantly, it is not trivial to get the correct effective index from the spectrum of an interfering device ^[13]. We need a quantitative discussion to choose the right effective index from many possible solutions.

In this paper, we address the above-mentioned challenges in a systematic way. As shown in Fig. 1, we first perform optical measurements on MZIs. We build a circuit model of the device and match the simulated transmission curve with measurement to get behavior parameters of a waveguide such as $n_{eff}$ and $n_{g}$ . Next, we build an accurate model to map $n_{eff}$ and $n_{g}$ to waveguide width and thickness. Then, using the model, we obtain waveguide geometry parameters from extracted behavior parameters. We automated the optical measurement and repeated the geometry extraction on devices with the same design over the wafer. From the extraction, we perform variability analysis and derive fabrication variation with high accuracy. In Section 2, we propose an improved geometry model to offer high modeling accuracy. In Section 3, we show that the curve fitting method is less sensitive to measurement noise and helps in removing GC envelope, and in Section 4, we show how to design two MZIs to extract the accurate and unique $n_{eff}$ and $n_{g}$ given the waveguide geometry variation range. Section 5 discusses the procedure of the waveguide geometry from a couple of MZIs. Finally, in Section 6, we apply the method to extract linewidths and thicknesses of SOI waveguides on a die fabricated by the IMEC multiproject wafer (MPW) service.

Work flow of extracting behavior parameters and fabricated geometry using optical measurements.

Figure 1.Work flow of extracting behavior parameters and fabricated geometry using optical measurements.

2. ACCURATE GEOMETRY MODEL

We cannot measure linewidth $w$ and thickness $t$ directly, so we infer them from the results of optical transmission measurements. In particular, we first extract the effective index $n_{eff}$ and group index $n_{g}$ of the waveguide. To calculate the fabricated waveguide linewidth $w$ and thickness $t$ from the effective index $n_{eff}$ and group index $n_{g}$ , we require a geometry model that links $n_{eff}$ and $n_{g}$ to $w$ and $t$ . To build the model, we simulated an oxide-clad Si waveguide cross section [Fig. 2(a)] with the COMSOL Multiphysics finite element method (FEM) solver. According to the IMEC technology handbook, a fabricated strip waveguide has a sidewall angle of 85°, so we simulated the waveguide with an 85° sidewall angle. The waveguide width is the bottom width of the trapezoid. We swept width from 440 nm to 500 nm and thickness from 195 nm to 235 nm, and calculated $n_{eff}$ and $n_{g}$ at 1550 nm wavelength. The linewidth-thickness grid in Fig. 2(b) can be mapped one-on-one to the simulated $n_{eff}$ – $n_{g}$ grid (black solid line with circles) in Fig. 2(c). Then, we fitted simulated $n_{eff}$ and $n_{g}$ using polynomials of $w$ and $t$ . We build first-, second-, and third-order polynomial models [dashed lines in Fig. 2(c)]. Both $n_{eff}$ and $n_{g}$ vary quite linearly with $w$ and $t$ . Nonetheless, the first-order fitted linear model shows a clear deviation (maximum 0.32% error in $n_{eff}$ and 0.40% in $n_{g}$ ) from simulated $n_{eff}$ and $n_{g}$ . The third-order model matches the simulation very accurately. Then, to get $w$ and $t$ from spectral measurements, we wrote $w$ and $t$ , as a polynomial of $n_{eff}$ and $n_{g}$ at 1550 nm as $w = p_{0} + Σ_{i = 1}^{n} p e_{i} n_{eff}^{i} + Σ_{j = 1}^{n} p g_{j} n_{g}^{j} + Σ_{j = 1, i = 1}^{m}, n p e g_{i j} n_{eff}^{i} n_{g}^{j},$ (1) $t = q_{0} + Σ_{i = 1}^{n} q e_{i} n_{eff}^{i} + Σ_{j = 1}^{n} q g_{j} n_{g}^{j} + Σ_{j = 1, i = 1}^{m}, n q e g_{i j} n_{eff}^{i} n_{g}^{j},$ (2)where $p_{0}$ and $q_{0}$ are constant terms, and $p e_{i}$ and $p g_{i}$ are coefficients in polynomials.

Figure 2.(a) Cross-section schematic of an oxide-clad SOI strip waveguide with a 85° sidewall angle; (b) width and thickness grid of strip waveguides; (c) effective and group indices of strip waveguides on the geometry grid using the COMSOL FEM simulation, and the first-, the second-, or the third-order polynomial mapping model.

Using a simulated model can introduce a simulation error, coming from a difference between the actual waveguide geometry (both dimensions, shape, and material properties) and the trapezoidal geometry model we used in the mode solver. We have considered the sidewall angle in our model, but still we do not know the actual geometry, so it is very hard to compensate for this error. This means that some parameters will be only relative. In addition, the mapping error is the difference between the simulated trapezoidal waveguide geometry and extracted trapezoidal waveguide geometry using the geometry model. To calculate the mapping error, we first use simulated $n_{eff}$ and $n_{g}$ for simulated geometry. Then, we calculate the waveguide geometry using a polynomial model with the simulated $n_{eff}$ and $n_{g}$ . The mapping error reduces significantly with the order of the polynomial (Table 1). The first-order model has a maximum error of several nanometers, which is comparable to the reported intra-wafer manufacturing variations in width of 0.78 to 2.65 nm, and in thickness of 0.83 to 4.16 nm ^[14–19]. Obviously, the modeling error is too large to extract intra-die variations (variations between the same devices on one die). For a good estimation of the fabricated geometry variation, especially to study variability on the intra-wafer level and the intra-die level, a lower modeling error is required. The third-order polynomial model has a maximum error of 0.05 nm for both linewidth and thickness (Fig. 3), which is one magnitude smaller than the fabrication variation. A low modeling error makes geometry extraction more accurate and credible.

	Error Δw [nm]	Error Δt [nm]
First order	5.10	4.57
Second order	0.14	0.18
Third order	0.06	0.08

Table 1. Error of Polynomial Models

View all Tables

Figure 3.Error contour plot of the proposed third-order polynomial model where $w$ ranges from 440 to 500 nm and thickness ranges from 195 to 235 nm. Left, width extraction error; right, thickness extraction error.

3. EXTRACTING neff AND ng FROM AN MZI USING THE CURVE FITTING METHOD

Before we can determine the waveguide geometry from $n_{eff}$ and $n_{g}$ , we need to measure those quantities on a device and isolate the values for a straight waveguide. This effectively means that we need a device where a transmission measurement can give us an accurate extraction of both $n_{eff}$ and $n_{g}$ of a straight waveguide only, and which is in the presence of measurement noise and variation in the coupling structures.

In Ref. ^[5], ring resonators were used to extract $n_{eff}$ and $n_{g}$ . A ring resonator has sharp resonance peaks that are easy to detect and less prone to error in peak positions. However, a ring uses either a bent waveguide or it combines bends with straight sections, and the round trip path also includes the coupling sections, which will also have different optical properties. As such, we cannot isolate $n_{eff}$ and $n_{g}$ of the straight or bent waveguide. The alternative is to use an MZI with two arms that are identical except for that one arm has a longer straight part than another ^[11]. Ideally, the $n_{eff}$ and $n_{g}$ of the two arms are identical so that the spectral response of the MZI is only dependent on the path length difference between two arms. We can use MZI to measure varying waveguide geometry under process variation at different locations. In practice, the path length difference in a single MZI is also induced by a difference in $n_{eff}$ in the two arms because of process variability. Also, there is also the difference in $n_{eff}$ that the bends will contribute to. Those differences also lead to the extraction error in $n_{eff}$ and $n_{g}$ . However, since the distance between waveguides within an MZI compact is within 100 μm, we can safely assume that the error will be much smaller than the device-to-device variation.

An interfering structure such as a ring or an MZI would have a constructive interference when interference order $m$ is an integer and $m = \frac{n_{eff, 0} Δ L}{λ_{res}},$ (3) $FSR = \frac{λ_{res}^{2}}{n_{g} Δ L},$ (4)where $n_{eff, 0}$ is the prior estimate at the resonance wavelength $λ_{res}$ and $Δ L$ is the physical path length difference. If we know the resonance order $m$ , we could get $n_{eff}$ by locating $λ_{res}$ in the output spectrum. It is natural to apply the resonance detection method ^[5,12] to locate resonances in the spectrum and get both $n_{eff}$ and $n_{g}$ . However, an MZI has a sinusoidal-like spectral transmission. Its curve is quite flat near a maximum or minimum. Especially when measurement noise is involved, it is hard to locate its resonance by the peak detection method (in Fig. 6, the green cross indicates detected valleys). Using only the maximum detection method leads to a significant error (Table 2) in effective index and group index extraction, making it not suitable for geometry extraction.

	neff	ng	Width [nm]	Thickness [nm]
Curve fitting using a GC model	2.319	4.291	466.0	211.8
Peak detection	2.318	4.302	462.0	213.8
Difference between two methods	0.001	0.009	4.0	2.0

Table 2. Comparison between the Peak Detection Method and the Curve Fitting Method^a

View all Tables

To improve the extraction accuracy, we used the curve fitting technique. It extracts parameters by minimizing the difference between a circuit simulation and the measurement data. While maximum/minimum extraction only uses information at the peaks and ignores information on the rest of the spectrum, the curve fitting method utilizes the information from the entire measured spectrum, which is more tolerant to the measurement noise and gives more reliable extraction.

We simulate the MZI circuit in Caphe ^[20], a circuit simulator that calculates the scattering matrix of the circuit from the defined scattering matrix of each component. We built a Caphe circuit model of the MZI in the same way as our device under test (DUT), which has two types of components: waveguide and multimode interferometer (MMI).

A. Eliminating the Effect of the Grating Couplers

For automated measurement, light is vertically coupled to the DUT using a pair of GCs. We need to remove the envelope of the GC before fitting the spectrum with the circuit model. We can remove the GC in two ways. In one way, we measure a pair of reference GCs, preferably close to the DUT. By subtracting the measured reference GC from the DUT spectrum in a log scale, we can normalize the transmission spectrum of the DUT. It is a very common practice in optical transmission measurements. However, this method is error-prone in the presence of fabrication variation in GCs, because the reference GC can be subtly different from the GCs connected to the DUT. Even more, the input and output GCs of the DUT can be different from each other. As shown in Fig. 5, MZI transmission after subtracting the reference is not we would expect: the linear-scale spectrum should be a sine-like curve with maxima of the same amplitude. The significant mismatch between the measured and fitted spectrum introduces a large error in the extracted parameters.

Figure 4.(a) Layout of the MZI under test. (b) Circuit schematic of the MZI.

Figure 5.We removed the GC envelope using a reference GC near the DUT. Fabrication variation caused the measured spectrum after GC removal to be far from ideal (as shown by the spectrum simulated by the circuit model), as ideally the peaks in the spectrum should have the same amplitude. After GC removal, we fitted the measured spectrum with the circuit model (Fig. 4), not including the GC. Red solid curve, measured transmission spectrum after removing the GC envelope using a reference GC. Blue dotted curve, fitted spectrum using the circuit model. Left, the low-order MZI. Right, the high-order MZI.

Figure 6.This figure shows the measured transmission spectrum (red solid curve) and fitted spectrum (blue dotted curve) using the circuit model including the polynomial GC model. Also, valleys of the spectrum (green cross) are found by the peak detection method. Left, the low-order MZI. Right, the high-order MZI.

Instead, we use a fourth-order polynomial to represent the logarithmic transmission spectrum of the combined input and output GC. We include two GCs in the circuit model together with the MZI. As shown in Fig. 6, the fitting is considerably enhanced where the simulation matches measurement nicely. From fitting, we can get circuit parameters such as effective index, group index, and coefficients of the polynomial describing the GC.

B. Fitting Accuracy versus MZI Order

In the transmission spectrum of the MZI, the positions of the peaks and valleys give information about the effective index $n_{eff}$ . The periodicity of the transmission spectrum is determined by the group index $n_{g}$ . An MZI with a low-order $m$ has only a few peaks/valleys in the measurement band, and therefore it will have a low accuracy on the extraction of $n_{g}$ . On the other hand, a high-order MZI can give a high accuracy of $n_{g}$ extraction. As explained by Dwivedi et al. in Ref. ^[11], the combination of a low-order MZI and a high-order MZI can give a good accuracy on both $n_{eff}$ and $n_{g}$ . When noise is mingled in the measured spectrum, it will induce in curve fitting an error. Then, it is a question of how fitting accuracy is related to the order of the MZI. In this work, we used the nonlinear least squares method to fit the transmission curves with a waveguide compact model. We built a circuit model of a waveguide with 470 nm width and 215 nm thickness and simulated the transmission spectrum. Then, we add a $\pm 0.2 dBm$ to the spectrum to emulate the typical measurement noise. Finally, we get $n_{eff}$ and $n_{g}$ using the curve fitting. The fitting error we presented is the estimate of 1.96 times standard deviations of each of the parameters, which provides confidence limits of approximately 95%. As shown in Table 3, we did not get an accurate $n_{g}$ from the 15-order MZI. This leads to a huge error in extracted width and thickness. The fitting error decreases with the interference order. Therefore, an increasing interference order of the MZI improves fitting accuracy.

Order	neff	ng	w [nm]	t [nm]	Errorneff	Errorng	Errorw	Errort
15	2.336264	4.280000	472.6998	212.7182	0.014705	0.718928	270.355700	137.651600
50	2.340747	4.290649	469.3341	215.6051	2.18×10−5	0.001366	0.511691	0.270922
100	2.339561	4.288268	470.0631	214.9168	1.07×10−5	0.000663	0.248367	0.130336
150	2.339729	4.288448	470.0189	214.9829	7.39×10−6	0.000460	0.172499	0.090586

Table 3. Fitting Error versus Interference Order

View all Tables

	Width w [nm]		Thickness t [nm]
Die Number	(X=0, Y=0)	(X=−2, Y=2)	(X=0, Y=0)	(X=−2, Y=2)
Mean, μ	469.34	463.88	211.58	212.86
Standard deviation, σ	1.26	1.18	0.30	0.42
Extraction error	0.37		0.26

Table 4. Statistical Results for the Manufacturing Variations of a 200 mm Wafer Fabricated through a 193 nm DUV Lithography Process

View all Tables

4. HOW TO DESIGN THE MZI FOR GEOMETRY EXTRACTION

When fitting the transmission curve of the MZI, the fitting algorithm implicitly assumes that the order $m$ of the MZI is sufficiently accurate, i.e., that the peak near the center wavelength of 1550 nm corresponds with the designed order $m$ . However, in the presence of fabrication variation, this is not necessarily the case, and as the designed order of the MZI increases, the uncertainty on the measured order increases. Therefore, the design parameters should be chosen such that the low-order MZI can be used to pin the order $m$ of the device unambiguously ^[11], and make a good estimate of $n_{eff}$ . The order of the high-order MZI should be chosen such that a maximum of information can be extracted, based on the estimate of $n_{eff}$ obtained from the low-order MZI. We discuss the design process for these devices.

A. Deciding the Interference Order under a Given neff Variation

Equation (4) shows that if we know the resonance order $m$ , we can calculate $n_{eff}$ from the peak locating $λ_{res}$ in the output spectrum. However, if fabrication variations can shift the spectrum more than half a free spectral range (FSR), we can no longer be certain of the order $m$ . Therefore, we should design the MZI with sufficiently low-order $m$ such that the order at the center wavelength is always within $m \pm 0.5$ , which means $\frac{n_{eff, 0} Δ L}{λ} = m, \frac{(n_{eff, 0} - Δ n_{eff} / 2) Δ L}{λ} > m - 0.5, \frac{(n_{eff, 0} + Δ n_{eff} / 2) Δ L}{λ} < m + 0.5 .$ Given the variation in $n_{eff}$ is $Δ n_{eff}$ , we can decide the order $m$ that fulfills the condition and would give us sufficient confidence: $m < \frac{n_{eff, 0}}{Δ n_{eff}} .$ (5)This is equivalent to a constraint on the length difference $Δ L$ between the two arms of the MZI: $Δ L < \frac{λ}{Δ n_{eff}} .$ (6)

B. Defining the Bounds for neff and ng from the Geometry Variation

Within the measurement interval, the spectrum of an MZI looks quasi-identical if we shift the interference order $m$ by an integer number. Without a proper confidence interval on $n_{eff}$ , there would be multiple solutions of $n_{eff}$ to fit the spectrum. As $n_{eff}$ and $n_{g}$ can be mapped to linewidth $w$ and thickness $t$ , we can derive the bound of $(n_{eff}, n_{g})$ from the confident interval of the geometry parameters $(w, t)$ , which are supplied by the fab. As presented in Section 2, $n_{g}$ and $n_{eff}$ can be accurately mapped on $w$ and $t$ by a third-order polynomial model. For simplicity of analysis in the derivation below, we use the linear geometry model, where $n_{eff} = n_{{eff}_{0}} + \frac{\partial n_{eff}}{\partial w} (w - w_{0}) + \frac{\partial n_{eff}}{\partial t} (t - t_{0}), n_{g} = n_{g_{0}} + \frac{\partial n_{g}}{\partial w} (w - w_{0}) + \frac{\partial n_{g}}{\partial t} (t - t_{0}) .$ When the bounds for linewidth and thickness form the rectangle $A B C D$ in $(w, t)$ space [Fig. 7(a)], the parameter range in $(n_{eff}, n_{g})$ space lies in between $n_{eff} \in [n_{{eff}_{1}}, n_{{eff}_{2}}]$ and $n_{g} \in [n_{g_{1}}, n_{g_{2}}]$ : $n_{{eff}_{1}} = n_{{eff}_{0}} + \frac{\partial n_{eff}}{\partial w} (w_{1} - w_{0}) + \frac{\partial n_{eff}}{\partial t} (t_{1} - t_{0}),$ (7) $n_{{eff}_{2}} = n_{{eff}_{0}} + \frac{\partial n_{eff}}{\partial w} (w_{2} - w_{0}) + \frac{\partial n_{eff}}{\partial t} (t_{2} - t_{0}),$ (8) $n_{g_{1}} = n_{g_{0}} + \frac{\partial n_{g}}{\partial w} (w_{2} - w_{0}) + \frac{\partial n_{g}}{\partial t} (t_{1} - t_{0}),$ (9) $n_{g_{2}} = n_{g_{0}} + \frac{\partial n_{g}}{\partial w} (w_{1} - w_{0}) + \frac{\partial n_{g}}{\partial t} (t_{2} - t_{0}),$ (10)where $w_{0}$ and $t_{0}$ are the nominal values for $w$ and $t$ variations. The range is a rectangle in $(w, t)$ space whose center is $(w_{0}, t_{0})$ . We assume $Δ w = w_{2} - w_{1}$ and $Δ t = t_{2} - t_{1}$ . As $(w, t)$ and $(n_{eff}, n_{g})$ follow a near-linear mapping, the vertices of the bound rectangle in the $(w, t)$ space $A$ , $B$ , $C$ , $D$ can be mapped to the $(n_{eff}, n_{g})$ space as $A^{'}$ , $B^{'}$ , $C^{'}$ , $D^{'}$ . Intuitively, the bound rectangle $A B C D$ in $(w, t)$ space is linearly transformed into a parallelogram $A^{'} B^{'} C^{'} D^{'}$ in $(n_{eff}, n_{g})$ space. The tilted boundary $A^{'} B^{'} C^{'} D^{'}$ is within the original rectangular boundary and much smaller. For the fundamental TE mode of an SOI oxide-clad waveguide, $\frac{\partial n_{g}}{\partial w}$ is negative while $\frac{\partial n_{eff}}{\partial w}$ , $\frac{\partial n_{eff}}{\partial t}$ , and $\frac{\partial n_{g}}{\partial t}$ are positive. The parallelogram $A^{'} B^{'} C^{'} D^{'}$ then will be tilted as in Fig. 7(b).

Figure 7.Bounds of the extraction. (a) The bound of width and thickness. (b) Rectangle bound [11] parallelogram, reduced bounds by linear transformation of geometry bounds. (c) Rectangle bounds cannot separate three groups of solutions (red, blue, and green circles). The parallelogram cleanly isolates the correct solutions (blue circles).

The bounds of the confidence interval for $n_{eff}$ are fundamental to get a correct extraction. Extraction of the group index $n_{g}$ does not pose that much of a problem, as the confidence interval is much larger, and there are no multiple solutions. Depending on whether we know $n_{g}$ of the same waveguide, we can estimate the range for $n_{eff}$ in two ways.

1. Estimating neff without Information on ng

Sometimes we cannot obtain accurate information of $n_{g}$ , such as when we are using a low-order MZI. Then, as in Ref. ^[11], we can calculate the uncertainty $Δ n_{eff}$ from geometry variations $Δ w$ and $Δ t$ as $Δ n_{eff, rectangle} = n_{{eff}_{2}} - n_{{eff}_{1}} = \frac{\partial n_{eff}}{\partial w} Δ w + \frac{\partial n_{eff}}{\partial t} Δ t,$ (11)which essentially corresponds to the width of the rectangle $A B C D$ , or the horizontal distance $A^{'} C^{'}$ .

2. Estimating neff with Information on ng

The maximal range of the $n_{eff}$ for a given $n_{g}$ is $E^{'} F^{'}$ [Fig. 7(b)], which is the maximal distance between two edges of the parallelogram at the given $n_{g}$ . The distance is dependent on the shape of the parallelogram. When $A^{'}$ is higher than $C^{'}$ [ $n (g, A) > n (g, C)$ ], using Eq. (7), we can derive $\frac{\partial n_{g}}{\partial t} Δ t < - \frac{\partial n_{g}}{\partial w} Δ w$ . Then, $E^{'} F^{'}$ is the horizontal distance between lines $A^{'} B^{'}$ and $C^{'} D^{'}$ , and the range of $n_{eff}$ is determined by the range of linewidth $Δ w$ . When $A^{'}$ is lower than $C^{'}$ [ $n (g, A) < n (g, C)$ ], using Eq. (7), we can derive $\frac{\partial n_{g}}{\partial t} Δ t > - \frac{\partial n_{g}}{\partial w} Δ w$ . Then, $E^{'} F^{'}$ is the horizontal distance between lines $A^{'} C^{'}$ and $B^{'} D^{'}$ , and the range of $n_{eff}$ is determined by the range of thickness $Δ t$ : $Δ n_{eff, parallelogram} = (- \frac{\frac{\partial n_{eff}}{\partial w} \frac{\partial n_{g}}{\partial t}}{\frac{\partial n_{g}}{\partial w}} + \frac{\partial n_{eff}}{\partial t}) Δ t, when \frac{\partial n_{g}}{\partial t} Δ t < - \frac{\partial n_{g}}{\partial w} Δ w,$ (12) $Δ n_{eff, parallelogram} = (\frac{\partial n_{eff}}{\partial w} - \frac{\frac{\partial n_{g}}{\partial w} \frac{\partial n_{eff}}{\partial t}}{\frac{\partial n_{g}}{\partial t}}) Δ w, when \frac{\partial n_{g}}{\partial t} Δ t > - \frac{\partial n_{g}}{\partial w} Δ w .$ (13)For the same geometry variation, an estimate of $n_{g}$ reduces the uncertainty $Δ n_{eff}$ . Figure 7(c) shows that we can separate three groups of solutions with the bound $A^{'} B^{'} C^{'} D$ , which are all located in the previous rectangle bound. Also, $Δ n_{eff}$ becomes solely dependent on $Δ t$ in the linear approximation. So if we have a small $Δ t$ and a relatively large $Δ w$ , we can benefit more from knowing $n_{g}$ , because we can use a high-order MZI that improves the accuracy of extraction. When $\frac{\partial n_{g}}{\partial t} Δ t < - \frac{\partial n_{g}}{\partial w} Δ w$ , the ratio between $Δ n_{eff}$ calculated without and with $n_{g}$ is $\frac{Δ n_{eff, rectangle}}{Δ n_{eff, parallelogram}} = a \frac{Δ w}{Δ t} + b,$ where $a = \frac{\frac{\partial n_{eff}}{\partial w}}{- \frac{\partial n_{eff}}{\partial w} \frac{\partial n_{g}}{\partial t} / \frac{\partial n_{g}}{\partial w} + \frac{\partial n_{eff}}{\partial t}}, b = \frac{\frac{\partial n_{eff}}{\partial t}}{- \frac{\partial n_{eff}}{\partial w} \frac{\partial n_{g}}{\partial t} / \frac{\partial n_{g}}{\partial w} + \frac{\partial n_{eff}}{\partial t}} .$ Both $a$ and $b$ are positive, so that the ratio is increasing with $\frac{Δ w}{Δ t}$ . Intuitively, the smaller that $A D$ is compared to $A B$ in Fig. 7(a), the shorter the $E F$ is. A similar conclusion can be made when $\frac{\partial n_{g}}{\partial t} Δ t < - \frac{\partial n_{g}}{\partial w} Δ w$ , where the ratio is increasing with $\frac{Δ t}{Δ w}$ .

5. EXTRACTION PROCEDURE

A. Extracting neff and ng from Two MZIs

The total process variation (intra-die, die-to-die, and wafer-to-wafer) on an isolated waveguide on an SOI platform is large. The variation can be several tens of nanometers for both linewidth and thickness. As discussed in the previous section, to capture the large variation using an MZI, we should choose a sufficiently low-order $m$ . However, this low-order MZI suffers from a low accuracy on $n_{g}$ extraction. On the other hand, a high-order MZI can offer good accuracy of $n_{g}$ extraction. So a combination of the two devices can give us both essential optical parameters. So we can extract $n_{eff}$ and $n_{g}$ using a low-order MZI and a high-order MZI: 1.

Extract a good estimate of $n_{eff}$ from the low-order MZI.

Extract an accurate $n_{g}$ from the high-order MZI.

Even though the devices are close together, they do not have the same $n_{eff}$ and $n_{g}$ because of local variations. To accurately map the waveguide geometry, we need to extract both $n_{eff}$ and $n_{g}$ from the same waveguide. In the following discussion, we will present how to obtain an accurate $n_{eff}$ and $n_{g}$ both of the high-order MZI step by step.

B. Extracting Inter-Die and Intra-Die Variability in Three Steps

In a wafer-scale fabrication process, we can identify different levels of process variations. For each die, all the variations that originate at levels such as lot-to-lot, wafer-to-wafer, and die-to-die variations have the same impact on every device in the die. In this paper, we categorize these variations together as the inter-die variation. However, then we also get the intra-die variation that affects devices differently on the same die. We can further decompose the intra-die variation into location-dependent variation and local variation. The location-dependent variation is the variation depending on the location of the device on the die. It can be caused by the continuous variation of thickness, photoresist spinning effects or plasma distributions, and other equipment nonuniformity that affects the fabricated geometry varied spatially. On the other hand, the local variation we define here induces local disparities between devices placed close together (less than a few hundred microns apart). It includes intrinsic variability such as thickness fluctuations and width randomness caused by pattern density nonuniformity. The sum of the three variations gives us the process variation of a device: ${variation}_{total} = {variation}_{inter-die} + {variation}_{location-dependent} + {variation}_{local} .$ (14)The total process variation is considerably larger than the intra-die variability. The total linewidth and thickness variation of an isolated waveguide on an SOI platform can amount to tens of nanometers ^[15], while intra-die variation is typically only a few nanometers ^[5,15,16,21].

With our two MZIs, we address variations on the different levels in three steps. In the first step, we extract $n_{eff}$ from a low-order MZI. Since the $n_{g}$ extracted from the low-order MZI is very inaccurate, we estimate the range of $n_{eff}$ without the information of $n_{g}$ by substituting geometry variation by the total variation in Eq. (11): $Δ n_{eff, total} = \frac{\partial n_{eff}}{\partial w} Δ w_{total} + \frac{\partial n_{eff}}{\partial t} Δ t_{total} .$ Given the range, we derive a fairly accurate value of $n_{eff}$ from the low-order MZI. In the second step, we obtain the $n_{eff}$ map over the die by interpolation. The map offers the average $n_{eff}$ of the waveguide placed at each location where we can remove the local variation, and the inter-die variation and location-dependent variation together determine the average value. In the third step, we use interpolated $n_{eff, μ}$ at each location as a reference. Now, rather than the total variation, we can only deal with the much smaller local variation. Since we can accurately extract $n_{g}$ on the high-order MZI, the range for $n_{eff}$ under the local variation is estimated by substituting geometry variation using the local variation in Eq. (12): $Δ n_{eff, local} = (- \frac{\frac{\partial n_{eff}}{\partial w} \frac{\partial n_{g}}{\partial t}}{\frac{\partial n_{g}}{\partial w}} + \frac{\partial n_{eff}}{\partial t}) Δ t_{local} .$ Because our analysis decreases the bound of extraction, we can use a much high-order MZI to get $n_{eff}$ and $n_{g}$ simultaneously and accurately.

C. Specification of the Two MZIs

The MZIs each consists of two waveguide arms and two 50-50 MMIs (Fig. 8 left). Our devices are fabricated by the IMEC MPW service. We design the waveguide with a linewidth of 450 nm. According to the technology handbook, the fabricated waveguide has a sidewall angle of 85°. The fabricated 450 nm waveguide will have a 470 nm mean value and $\pm 20 nm$ variations. The nominal thickness is 215 nm, and the variation is $\pm 10 nm$ . The pre-estimated $n_{eff, 0}$ of the $470 nm \times 215 nm$ waveguide is 2.340. Using Eq. (11), we calculated that the variation of $n_{eff, total}$ is $Δ n_{eff, total} = \frac{\partial n_{eff}}{\partial w} Δ w_{total} + \frac{\partial n_{eff}}{\partial t} Δ t_{total} = 0.002055 \times 40 nm + 0.003916 \times 20 nm = 0.16 .$ From the third-order model, we calculate that the $Δ n_{eff, total}$ is also 0.16. The arm length difference of the low-order MZI is $Δ L < \frac{1.55 μm}{0.16} = 9.7 μm$ , and the order of $m_{l} < \frac{n_{eff, 0}}{Δ n_{eff, total}} \approx 15.00$ .

Figure 8.Top left, low-order and high-order MZIs we used for geometry extraction. Bottom left, locations of two devices on a die. Right, locations of dies on the wafer. Red grid indicates dies on the wafer. The black circle is the boundary of the wafer.

Figure 9.Extracted $n_{eff}$ and $n_{g}$ of the high-order MZI. Left, die ( $X = 0$ , $Y = 0$ ); right, die ( $X = - 2$ , $Y = 2$ ).

Only a few references report typical fabricated geometry maps of silicon waveguides on die-level. Thickness maps in SOI depend largely on the qualities of the source wafer. Linewidth maps depend much more on the actual fabrication process and will be very different for devices fabricated with deep UV lithography or e-beam lithography, and vary between fabs. From the research of Lu et al. ^[5], the thickness varies slowly over the die, and the maximum difference between neighboring thickness is 0.5 nm. So we assume as a worst case that the thickness is slowly changing locally with only $Δ t_{local} = 2.00 nm$ . We also assume that $Δ w_{local} = 15.00 nm$ , which is significant. Using Eq. (12), we calculated the variation of $n_{eff, intra-die}$ . Notice that we have no pre-estimate of the nominal value of $w$ and $t$ locally so that they can be any value under the total variation. We calculated $- \frac{\frac{\partial n_{eff}}{\partial w} \frac{\partial n_{g}}{\partial t}}{\frac{\partial n_{g}}{\partial w}} + \frac{\partial n_{eff}}{\partial t}$ for $w \in [450, 490] nm$ and $t \in [205, 225] nm$ . Its value is in between 0.0045 and 0.0074. Therefore, the max local $n_{eff}$ variation we can surely extract is $Δ n_{eff, intra-die} = (- \frac{\frac{\partial n_{eff}}{\partial w} \frac{\partial n_{g}}{\partial t}}{\frac{\partial n_{g}}{\partial w}} + \frac{\partial n_{eff}}{\partial t}) Δ t_{intra-die} = 0.0074 \times 2 = 0.0148,$ which is the same calculated using the third-order model. The arm length difference of the high-order MZI is $Δ L_{h} < \frac{λ}{Δ n_{eff, local}} = 104.7 μm$ , where the order is around $m \approx 158$ . The confidence limit estimated with the information on $n_{g}$ on the order is $2.93 \times$ that calculated by the original method without the information on $n_{g}$ . Based on the estimation, we design the low-order MZI to have an order around $m = 15$ at 1550 nm and the order of the high-order MZI is $m = 150$ .

6. RESULTS

We automated optical measurements on 21 dies in the same wafer. The optical measurement was conducted in our clean room with the room temperature controlled at 20°C. On each die, we distributed 44 copies of the MZI pair (Fig. 8, bottom left) and repeated the fitting for all MZI blocks. With the estimated process variations on different levels, we set up the bound for $n_{eff}$ and followed the three-step procedure to extract $n_{eff}$ and $n_{g}$ of high-order MZIs (Fig. 6, right). Each point in the scatter plot represents the extracted $n_{eff}$ and $n_{g}$ of one waveguide on the die (Fig. 9). All the points are gathered in one group as confined by the bound. The average fitting error for $n_{eff}$ is $1.1 \times 10^{- 5}$ , and the average fitting error for $n_{g}$ is $1.0 \times 10^{- 3}$ . These fitting errors propagate to fitting errors of 0.33 nm in width and 0.18 nm in thickness. Adding the mapping error of the geometry model in Section 2, the total extraction errors for width $w$ and thickness $t$ are ${Error}_{w} = 0.31 + 0.06 = 0.37 nm, {Error}_{t} = 0.18 + 0.08 = 0.26 nm .$ Using the geometry model, we mapped $n_{eff}$ and $n_{g}$ to width $w$ and thickness $t$ of the high-order MZI arms. The extracted linewidth on the die ( $X = 0$ , $Y = 0$ ) in the wafer center (Fig. 8, right) ranges from 468.8 nm to 471.9 nm, and the thickness ranges from 211.4 nm to 212.3 nm. The standard deviations are 1.26 nm and 0.30 nm, respectively (Table 4). The extracted linewidth on the die ( $X = - 2$ , $Y = 2$ ) near the boundary of the wafer ranges from 461.4 nm to 466.8 nm, and the thickness ranges from 212.1 nm to 214.0 nm. The standard deviations are 1.18 nm and 0.42 nm, respectively (Table 4). For both dies, we observed a very weak correlation (correlation coefficient is $- 0.2856$ ) between the linewidth and the thickness. We fitted width and thickness to its location on the die with a linear model, and the green grid is the fitted map (Fig. 10) that indicates the location-dependent variation. We did observe a systematic trend for width $w$ on this die, but the trend is quite flat on each die, which is less than 1 nm. Meanwhile, the width $w$ shows an obvious local variation with a maximum of 3 nm. For thickness $t$ , both dies exhibit location dependency, which might be the result of slow-varying systematic variation over the wafer. Local variation has a maximum of 0.5 nm, well below the local variation range we set in the extraction. The maximum location-dependent difference in thickness on the die ( $X = 0$ , $Y = 0$ ) is 0.4 nm while on the die ( $X = - 2$ , $Y = 2$ ) is 1.5 nm.

Figure 10. $x$ and $y$ coordinates give the locations of the MZIs on two dies. Blue solid dot, extracted value. Green grid, fitted map of extracted values using a bivariate polynomial. (a) Extracted width map of die ( $X = 0$ , $Y = 0$ ) (in the center of the wafer). (b) Extracted thickness map of die ( $X = 0$ , $Y = 0$ ). (c) Extracted width map of die ( $X = - 2$ , $Y = 2$ ) (near the edge of the wafer). (d) Extracted thickness map of die ( $X = - 2$ , $Y = 2$ ).

We also extracted a width $w$ and thickness $t$ wafer map of one pair of MZIs that share the same location on every die. The extracted wafer map (Fig. 11) shows an explicit location dependence of fabricated geometry. We fit a parameter wafer map of $w$ , $t$ using a second-order bivariate polynomial. The slow-varying trend of the linewidth matches the dome-like radial symmetric pattern of the wafer-level systematic variation. The width systematic variation ranges from 459 nm to 465 nm while the random part has a maximum 2 nm contribution. The thickness also presents a strong location dependence. Its systematic variation ranges from 211.5 nm to 214.5 nm. Its random variation has a maximum of less than 1 nm.

Figure 11.We extracted the linewidth and thickness on the same device over 21 dies on the wafer. Top left, systematic linewidth variation; bottom left, random linewidth variation; top right, systematic thickness variation; bottom right, random thickness variation.

7. CONCLUSION

We show how to extract waveguide geometry from optical transmission measurement, and discuss how to increase the accuracy compared to existing methods. We replaced the linear mapping model between $(w, t)$ and $(n_{eff}, n_{g})$ with an accurate third-order geometry model to obtain accurate waveguide geometry from its effective index and group index. We applied a curve fitting method that is less sensitive to measurement noise and helps in removing the GC envelope. We discussed how to set parameter bounds under a given process variation, which helps to choose the correct set of extracted parameter values from multiple solutions. We showed that, when the group index of a waveguide is known, we can reduce the parameter bounds for the effective index, allowing us to use a higher order MZI and improve fitting accuracy. We proposed a procedure to separate different levels of process variation so that our method can deal with a total variation of several tens of nanometers and still obtain accurate linewidth and thickness extraction. We applied the method to measurement data from two dies and presented the linewidth and thickness map on die-level. We also applied the method to extract one pair of MZIs in 21 dies and presented a simple wafer map of fabricated geometry.

Acknowledgment

Acknowledgment. Part of this work is supported by the Flemish Research Foundation (FWO-Vlaanderen) and by the Flemish Agency for Innovation and Entrepreneurship (VLAIO) with the MEPIC project.

APPENDIX A: PROCEDURE TO APPLY THE CURVE FITTING M

Calibrate the laser before the measurement to assure the wavelength accuracy of the measurement spectrum.

Calibrate the laser after the entire measurement to assure the stability of the laser.

Measure a reference pair of GCs. Fit its log-scale transmission with a four-order polynomial. Use the fitted coefficients of the polynomial as the nominal values for GC parameters in the circuit model.

Use the circuit model of the MZI and the curve fitting method to extract $n_{eff}$ from the low-order MZI. Estimate the bound of $n_{eff}$ from the total fabrication variation on the waveguide linewidth and thickness. Since we have no accurate knowledge of $n_{g}$ of the low-order MZI, set the bound of $n_{eff}$ using Section 4.B.1.

Extract $n_{eff}$ of the low-order MZI over the die. Fit the $n_{eff}$ die map using first-order/second-order polynomials.

Use the circuit model of the MZI and the curve fitting method to extract $n_{eff}$ and $n_{g}$ from the high-order MZI. Estimate the nominal value of $n_{eff}$ at each die location using the fitted $n_{eff}$ die map from low-order MZIs. Estimate the bound of $n_{eff}$ from intra-die fabrication variation on the waveguide linewidth and thickness. Since we can extract $n_{g}$ accurately from the high-order MZI, set the bound of $n_{eff}$ using Section 4.B.2.

Extract $n_{eff}$ and $n_{g}$ of the high-order MZI over the die.

Map $n_{eff}$ and $n_{g}$ to linewidth $w$ and thickness $t$ using the third-order geometry model.

Plot the die map of width and thickness. Identify obvious outliers in the thickness map and check the errors of the corresponding curve fitting to evaluate the problem.

References

[1] S. Pathak, D. Van Thourhout, W. Bogaerts. Design trade-offs for silicon-on-insulator-based AWGs for (de)multiplexer applications. Opt. Lett., 38, 2961-2964(2013).

[2] A. Ribeiro, S. Dwivedi, W. Bogaerts. A thermally tunable but athermal silicon MZI filter. 18th European Conference on Integrated Optics 2016 (ECIO)(2016).

[3] S. Dwivedi, H. D’heer, W. Bogaerts. Maximizing fabrication and thermal tolerances of all-silicon FIR wavelength filters. IEEE Photon. Technol. Lett., 27, 871-874(2015).

[4] W. A. Zortman, D. C. Trotter, M. R. Watts. Silicon photonics manufacturing. Opt. Express, 18, 23598-23607(2010).

[5] Z. Lu, J. Jhoja, J. Klein, X. Wang, A. Liu, J. Flueckiger, J. Pond, L. Chrostowski. Performance prediction for silicon photonics integrated circuits with layout-dependent correlated manufacturing variability. Opt. Express, 25, 9712-9733(2017).

[6] T.-W. Weng, Z. Zhang, Z. Su, Y. Marzouk, A. Melloni, L. Daniel. Uncertainty quantification of silicon photonic devices with correlated and non-Gaussian random parameters. Opt. Express, 23, 4242-4254(2015).

[7] Y. Xing, D. Spina, A. Li, T. Dhaene, W. Bogaerts. Stochastic collocation for device-level variability analysis in integrated photonics. Photon. Res., 4, 93-100(2016).

[8] W. Bogaerts, M. Fiers, P. Dumon. Design challenges in silicon photonics. IEEE J. Sel. Top. Quantum Electron., 20, 8202008(2014).

[9] N. Ayotte, A. D. Simard, S. Larochelle. Long integrated Bragg gratings for SOI wafer metrology. IEEE Photon. Technol. Lett., 27, 755-758(2015).

[10] X. Wang, W. Shi, H. Yun, S. Grist, N. A. F. Jaeger, L. Chrostowski. Narrow-band waveguide Bragg gratings on SOI wafers with CMOS-compatible fabrication process. Opt. Express, 20, 15547-15558(2012).

[11] S. Dwivedi, A. Ruocco, M. Vanslembrouck, T. Spuesens, P. Bienstman, P. Dumon, T. Van Vaerenbergh, W. Bogaerts. Experimental extraction of effective refractive index and thermo-optic coefficients of silicon-on-insulator waveguides using interferometers. J. Lightwave Technol., 33, 4471-4477(2015).

[12] T. Horikawa, D. Shimura, H. Takahashi, J. Ushida, Y. Sobu, A. Shiina, M. Tokushima, S.-H. Jeong, K. Kinoshita, T. Mogami. Extraction of SOI thickness deviation based on resonant wavelength analysis for silicon photonics devices. IEEE SOI-3D-Subthreshold Microelectronics Technology Unified Conference (S3S), 10, 1-3(2017).

[13] W. Bogaerts, L. Chrostowski. Silicon photonics circuit design: methods, tools and challenges. Laser Photon. Rev., 12, 1700237(2018).

[14] S. K. Selvaraja, W. Bogaerts, P. Dumon, D. Van Thourhout, R. Baets. Subnanometer linewidth uniformity in silicon nanophotonic waveguide devices using CMOS fabrication technology. IEEE J. Sel. Top. Quantum Electron., 16, 316-324(2010).

[15] S. K. Selvaraja, E. Rosseel, L. Fernandez, M. Tabat, W. Bogaerts, J. Hautala, P. Absil. SOI thickness uniformity improvement using corrective etching for silicon nano-photonic device. IEEE International Conference on Group IV Photonics GFP, 71-73(2011).

[16] S. K. Selvaraja. Wafer-Scale Fabrication Technology for Silicon Photonic Integrated Circuits(2011).

[17] R. G. Beausoleil, A. Faraon, D. Fattal, M. Fiorentino, Z. Peng, C. Santori. Devices and architectures for large-scale integrated silicon photonics circuits. Proc. SPIE, 7942, 794204(2011).

[18] X. Chen, M. Mohamed, Z. Li, L. Shang, A. R. Mickelson. Process variation in silicon photonic devices. Appl. Opt., 52, 7638-7647(2013).

[19] 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, P. Absil. 193??nm immersion lithography for high-performance silicon photonic circuits. Proc. SPIE, 9052, 90520F(2014).

[20] M. Fiers, T. Van Vaerenbergh, J. Dambre, P. Bienstman. CAPHE: time-domain and frequency-domain modeling of nonlinear optical components. Advanced Photonics Congress, IM2B.3(2012).

[21] D.-X. Xu, J. H. Schmid, G. T. Reed, G. Z. Mashanovich, D. J. Thomson, M. Nedeljkovic, X. Chen, D. Van Thourhout, S. Keyvaninia, S. K. Selvaraja. Silicon photonic integration platform—have we found the sweet spot?. IEEE J. Sel. Top. Quantum Electron., 20, 8100217(2014).