Modeling and analysis of sulfur hexafluoride plasma etching for silicon microcavity resonators

The microcavities are fabricated using a two-step plasma etching process which is compatible with CMOS technology [10]. The process is performed on a moderately doped, n-type (100) $\mathrm{Si}$ wafer. Three layers of AZ6624 photoresist are added to a total thickness of 9 µm. The photoresist is then lithographically patterned with an array of 100 cylindrical holes with linearly increasing diameters between 12.4 µm and 52 µm.

The processing then follows the procedure illustrated in figure 1(a)). The wafer is etched for 320 s in an inductively coupled plasma (ICP) of $\mathrm{SF}_6$ gas at a flow rate of 100  sccm, a substrate temperature of 30 ^∘C, an ICP power of 2 kW, and a table power of 15 W. The photoresist is afterward completely removed with acetone. The second etch step follows the same recipe and is performed for 48 min.

To reduce the surface roughness, wet oxidation is applied for smoothing. For this purpose, the wafer is subjected to wet oxidation to a thickness of $2~\mu$m, followed by oxide removal using a solution of hydrogen fluoride ($\mathrm{HF}$). This polishing procedure is performed twice.

From the set of manufactured microcavities, we select three representative samples. They have an initial photoresist opening diameter d of 12.4 µm, 34 µm, and 52 µm. Their topography is recorded after the polishing procedure with a white light profilometer.

Near the edges of the microcavities, the angles are steep, leading to reflections and therefore loss of profilometer topography data. In order to gather information even in these regions of the specimen, a second white light profilometer measurement is performed on the tilted wafer. The two measurements are processed and tilt-corrected using Gwyddion [33]. After removal of the lost points due to reflection and matching the centers of the cavities, the measurements are combined into a single list of coordinates for each of the three microcavities.

The optical characteristics of resonators are fundamentally defined by their geometry. This motivates the use of an accurate simulation tool to capture the time evolution of the surface during the whole fabrication process. This is achieved by employing the topography simulator ViennaTS [34]. This simulator uses the level-set method [35], which describes the evolving surface as the zero level-set of the signed distance function $\phi(\boldsymbol{x)}$. The time evolution of the surface is then described by the level-set equation [19]

$$\begin{equation} \frac{\partial \phi (\boldsymbol{x}, t)}{\partial t} + V(\boldsymbol{x})|\nabla \phi(\boldsymbol{x}, t)| = 0 \, , \end{equation} \tag{ 1 }$$

where $V(\boldsymbol{x})$ is the scalar velocity field describing the local etch or deposition rates.

The topography simulator uses a velocity field $V(\boldsymbol{x})$ defined by a pseudo-particle model to represent the flux of reactants to the surface [23, 36]: the complex plasma species produced in the reactor are abstracted into pseudo-particles which interact with the evolving surface. These particles are generated on a source plane above the evolving surface according to a specified particle distribution and subsequently tracked via a Monte Carlo ray-tracing method [37]. Upon interacting with the surface, the pseudo-particles contribute to the local velocity field and subsequently are re-emitted according to a prescribed model.

Despite the use of pseudo-particles, our model represents a continuum approximation of the materials and their processing. The pseudo-particles do not represent individual reactant species such as ions or molecules, instead, they are Monte Carlo samples of the continuum flux integral [37]. The model is therefore unable to reproduce surface roughness caused by local inhomogeneity on the molecular scale [38]. However, due to the intrinsic randomness of the pseudo-particles, some artificial roughness can be present from the Monte Carlo noise. Our simulation approach does not incorporate any smoothing procedure, instead, we rely on having a sufficiently large number of samples in order to keep noise under control.

Our topography simulation approach seeks to replicate each experimental processing step presented in section 2.1. The initial geometry is a 160 µm × 160 µm slab of $\mathrm{Si}$ covered with a photoresist layer with a thickness of 9 µm with a cylindrical opening, as presented in figure 1(b)). Both plasma etching steps are represented by the same pseudo-particle model described in section 3.2 with different etch rate parameters for each step. The mask removal step is modeled by geometrically removing the entire photoresist, that is, it assumes that the removal step is ideal. Since our model does not include surface roughness, the only consequence of the polishing procedure described in section 2.1 is additional etching of the whole wafer surface. Our model captures this behavior by artificially increasing the etch rate for the second plasma etch step.

At the end of the simulation (cf figure 1(c)), the level-set is converted to a triangle mesh (VTK format [39]). This mesh is then analyzed according to the procedure presented in section 4.

The pseudo-particle model, illustrated in figure 2, tracks individual Monte Carlo pseudo-particles generated with an isotropic distribution on a source plane above the evolving surface with reflective boundary conditions. All pseudo-particles are of the same type and have the same properties. A pseudo-particle interacts with both the photoresist and the $\mathrm{Si}$ by etching and subsequently reflecting with an isotropic re-emission distribution.

Zoom In Zoom Out Reset image size

This model is uniquely determined by designating a plane wafer etch rate (PWR) and sticking probability β for each material being either photoresist ($\mathrm{res}$) or $\mathrm{Si}$. In general, sticking probabilities are complex quantities which vary locally, depending on, e.g. the clean surface sticking probability and surface coverages by reactants [21, 32, 40, 41]. For simplicity, we approximate the sticking probability as a constant for each material. This approximation is well-established in the literature and has shown agreement with experiment [42]. Upon intersecting with the surface at point x , each particle p adds to the local scalar velocity $V_\mathrm p(\boldsymbol{x})$ according to

$$\begin{equation} V_\mathrm p(\boldsymbol{x}) = -p_\mathrm{p} \times F_{p} \times \mathrm{PWR}_{\mathrm{Si}/\mathrm{res}} \, , \end{equation} \tag{ 2 }$$

where $F_\mathrm{p}$ is the per-particle flux and $p_\mathrm{p}$ a per-particle probability. Since p is a Monte Carlo pseudo-particle, its corresponding $p_\mathrm p$ can be interpreted as a weight to the pseudo-particle payload $F_\mathrm{p} \times \mathrm{PWR}_{\mathrm{Si}/\mathrm{res}}$. After each reflection, said weight is reduced according to the fixed sticking coefficient. Thus, $p_\mathrm p$ is initialized as 1 for each pseudo-particle generated in the source plane. Upon encountering a surface element, p is terminated and a new, reflected pseudo-particle is generated with probability

$$\begin{equation} p_{\mathrm{refl}} = \beta_{\mathrm{Si}/\mathrm{res}} \times p_\mathrm{p}. \end{equation} \tag{ 3 }$$

The final etch rate is then obtained by combining all particle contributions to the velocity field:

$$\begin{equation} V(\boldsymbol{x}) = \sum_p V_\mathrm p(\boldsymbol{x}). \end{equation} \tag{ 4 }$$

As we show in the following, this model captures the essential characteristics of the low-voltage-bias $\mathrm{SF}_6$ etch process described in section 2.1. $\mathrm{SF}_6$ is known to be a source of highly reactive, non-selective F atoms [43, 44] which chemically etch the surface both vertically and laterally [17]. This motivates our use of an isotropic description of the source and re-emission distributions. It also highlights the importance of visibility and reflection effects. These effects underscore the difference in reactant transport between low-voltage-bias plasma etching and isotropic wet etching processes which are also common in micromachining processing [45]. Wet etching processes in the reaction-limited regime are expected to have an equal supply of reactants to the exposed geometry [46], ultimately leading to a different surface.

The pseudo-particle model described in section 3.2 has the plane wafer etch rates PWR$_{\mathrm{Si}/\mathrm{res}}$ and sticking probabilities $\beta_{\mathrm{Si}/\mathrm{res}}$ as free input parameters. Additionally, we allow the plane wafer etch rate for Si PWR$_{\mathrm{Si}}$ to be different for each etch step, enabling the capture of reactor loading effects. The model thus has five free parameters in total, which are described in table 1. Although some of these parameters can be experimentally estimated we assume that they are free to achieve maximum flexibility. That being said, having five free parameters, a robust calibration procedure to match the simulated geometries to the experimental data, as is illustrated in figure 3, is required.

Zoom In Zoom Out Reset image size

Table 1. Definition of the five free pseudo-particle model parameters.

Free parameter	Definition
1st PWR$_{\mathrm{Si}}$	Plane wafer etch rate of $\mathrm{Si}$ for 1st step
2nd PWR$_{\mathrm{Si}}$	Plane wafer etch rate of $\mathrm{Si}$ for 2nd step
PWR$_{\mathrm{res}}$	Plane wafer etch rate of the photoresist
$\beta_{\mathrm{Si}}$	Sticking probability of $\mathrm{Si}$
$\beta_{\mathrm{res}}$	Sticking probability of the photoresist

The central challenge of calibrating a topography simulation to experimental data is that of matching geometrical data. From the simulation (cf section 3.1 and figure 3(c)) we obtain a triangle mesh, whereas the characterization procedure described in section 2.2 provides a list of coordinates. We match such disparate data by developing a purpose-built feature detection algorithm and applying it to both experimental and simulated geometries, as shown in figure 3(d)). The underlying method is the circular Hough transform [47], as implemented in the OpenCV library [48]. This method enables the detection of the circle which forms the maximum opening of the experimental and simulated microcavities, delimiting the relevant area. All of the data points are then projected into their radial coordinate, reducing the data set from three to two dimensions. This projection is possible since rotational symmetry is present. We then obtain a functional description of the microcavities by fitting an even sixth-order polynomial to all points on the inside. As summarized in figure 3(e)), after the feature detection procedure, each simulated and experimental geometry is described by their maximum opening diameter O, maximum depth h, and by a mirror shape given by a polynomial $\mathrm{poly}(x)$.

The calibration procedure requires the definition of a cost function to be minimized, illustrated in figure 3(g)). We define our cost function as the Euclidean norm of a vector of residuals. For each microcavity to be calibrated, we construct the residuals as

$$\begin{equation} \left( \begin{array}{c} |h_{\mathrm{exp}} - h_{\mathrm{sim}}| \\ \frac{1}{N_x}\sum_x |\mathrm{poly}_{\mathrm{exp}}(x) - \mathrm{poly}_{\mathrm{sim}}(x)| \end{array} \right) \, , \end{equation} \tag{ 5 }$$

where N_x is the number of samples of x taken from the region inside the microcavity. We calibrate our simulation ($\mathrm{sim}$) simultaneously to m multiple experimentally ($\mathrm{exp}$) realized microcavities, each with a different photoresist opening d. Since each microcavity contributes with two entries to the vector of residuals, the final vector contains $2\,\times\,m$ entries.

In order to determine the global minimum of the cost function, the generalized simulated annealing (GSA) algorithm is employed [49], as implemented in the SciPy library [50, 51] (cf figure 3(h)). Since our topography model is based on a Monte Carlo method [37], our cost function is susceptible to random variations after slight changes in input parameters. That is, our minimization landscape is noisy with respect to local searches in the parameter space. Therefore, we employ the GSA algorithm to minimize our cost function. The GSA algorithm enables larger jumps in the parameter space and searches more efficiently, in our case, compared to other algorithms available in SciPy.

The calibration approach is performed, as discussed in section 4.1, on m = 3 microcavities. The lateral extent of the simulation domain is 160 µm × 160 µm with a finite difference grid spacing of 0.5 µm. The number of Monte Carlo pseudo-particles per time step is of approximately 6.1 million (120 pseudo-particles per source plane grid point). The optimization is performed using a dual-socket compute node equipped with dual Intel Xeon Gold 6248 processors (40 physical computing cores in total). Since each mirror can be calculated independently and simultaneously, each simulation is assigned 13 threads in order to evenly distribute the workload over the 40 computing cores. The total runtime depends greatly on the chosen parameter value bounds, particularly on the plane wafer etch rates. A typical GSA run with 100 iterations requires approximately 60 hours.

Following the procedure described in section 3.3, we calibrated our pseudo-particle model to three microcavities with d of: 12.4 µm, 34 µm, and 52 µm. A visual comparison of the calibrated simulation to the profilometer data is summarized in figure 4. The final geometries of each of the three microcavities, obtained by the profilometer according to the procedure described in section 2.2, are the calibration targets for the residuals constructed according to (5). The calibrated parameters are shown in table 2.

Zoom In Zoom Out Reset image size

Table 2. Calibrated parameters.

Free parameter	Calibrated value
1st PWR$_{\mathrm{Si}}$ (µm min−1)	2.151
2nd PWR$_{\mathrm{Si}}$ (µm min−1)	0.661
PWR$_{\mathrm{res}}$ (µm min−1)	0.209
$\beta_{\mathrm{Si}}$ (-)	$7.5\%$
$\beta_{\mathrm{res}}$ (-)	$6.1\%$
Cost function (µm)	0.367

The calibrated PWR$_{\mathrm{Si}/\mathrm{res}}$ and $\beta_{\mathrm{Si}/\mathrm{res}}$ shown in table 2 are consistent with values reported from the literature [45, 52–54]. The $70\%$ reduction of PWR$_\mathrm{Si}$ during the second etch step is consistent with the loading effect [17, 55, 56]. After the removal of the photoresist, there is a larger available area of the $\mathrm{Si}$ to be consumed by the reactants. This locally reduces the availability of reactants, decreasing the etch rate. Additionally, a photoresist selectivity of approximately 10 to 1 is consistent with reported results [52, 57, 58].

As described in section 1, the goal of manufacturing these cavities is the production of optical resonators [10]. For our analysis, we therefore consider that the resonator is assembled in a plano-concave geometry (cf figure 5). Two additional parameters are then extracted from the simulated geometry arising from the calibrated parameters presented in section 4.1: the ROC and the beam waist size $w_\mathrm M$. The extraction procedure follows the scheme presented in figure 5.

Zoom In Zoom Out Reset image size

Firstly, the feature detection algorithm presented in section 3.3 outputs a detected opening O. From O, we extract a representative ROC of the resonator by fitting a parabola $y = ax^2 + c$ in the region inside $50\%$ of the final opening. By focusing on a central and relatively small region, the extracted parameters are accurate and still representative of the global characteristics. The ROC is then:

$$\begin{equation} \mathrm{ROC} = \frac{1}{2|a|}. \end{equation} \tag{ 6 }$$

For resonators in the plano-concave configuration, the Gaussian beam waist size at the concave mirror $w_\mathrm M$ is given by [59]

$$\begin{equation} w_\mathrm M = \sqrt{\frac{\lambda}{\pi}\sqrt{\frac{L \times (\mathrm{ROC})}{1-\frac{L}{\mathrm{ROC}}}}} \, , \end{equation} \tag{ 7 }$$

where L is the resonator length and λ a representative wavelength. Here, we assume that the ratio $L/(\mathrm{ROC}) = 0.75$ and that the relevant wavelength is $\lambda = 1.55~\mu$m, which are the values obtained from further infrared laser analysis of assembled resonators [10].

The $w_\mathrm M$ can be interpreted in a Gaussian sense as one standard deviation σ of the beam relative to its center of symmetry. This motivates us to define the region of interest as $6w_\mathrm M$, or equivalently, 3σ. This represents enough area such that the maximum possible finesse, which is a quantity to measure losses inside the resonator, is larger than 10⁷ [60, 61]. Therefore, we extract a functional description of the microcavity by fitting a sixth-order even polynomial to the region inside either of $6w_\mathrm M$ or the detected opening O, whichever is smallest.

Based on the calibrated parameter set different etch time scenarios can be accurately investigated. The devices fabricated as described in section 2.1 were etched with the photoresist for $t_{\mathrm{resist}} = 320$ s for a cumulative etch time $t_{\mathrm{etch}} = 3200$ s. Our approach allows for investigating possible states of the topography for different etch times, not only for intermediate periods during the executed etch but also for alternative etch time configurations. This highlights how our modeling approach is able to provide detailed insights which would have otherwise required numerous additional experiments. This approach is valid under the moderate assumptions that the reactor conditions, as reflected in the plane-wafer etch rates and sticking probabilities, are constant during each step [28].

The evolution of selected geometric features of a microcavity with a photoresist $d = 12.4~\mu$m is depicted in figure 6. We observe that $t_{\mathrm{resist}}$ directly defines the h of the microcavity. This is expected, as after the removal of the photoresist the whole wafer is exposed to the reactor and there is comparatively little variation in local etch rates. The opening O is influenced by both etching steps; however, the rate of increase of O is influenced by $t_{\mathrm{resist}}$. Interestingly, the ROC has its minimum defined by the first etch step, reaching a global minimum near $t_{\mathrm{resist}} = 1000$ s. Since a small ROC is a desired attribute of the device, this indicates that the second etch step should be kept as short as possible.

Zoom In Zoom Out Reset image size

The noise observed in the extracted geometrical features is a consequence of our simulation engine using a Monte Carlo method [37]. Nonetheless, we highlight that our simulations match the experimental data when applying the same geometrical feature extraction procedure described in section 4.1 to the experimental topography.

Scale analysis is a procedure to determine whether the process evolution characteristics scale with a particular geometric parameter. We are interested in determining if the processing characteristics scale with the main geometric input: the photoresist d. We achieve this by dividing both the etch time $t_{\mathrm{resist}}$ and an output variable, such as the ROC, by d. With the scale analysis, we can evaluate if higher or lower values of d lead to resonators of higher quality.

As it can be seen in figure 7, our simulations show that the scaling curves for all three microcavities have similar shapes. This shows that the process, under the conditions described in section 2.1, scales directly with d. Therefore, we expect that other geometrical scales, such as the thickness of the photoresist or the aspect ratio of its cylindrical hole, do not play a significant role, otherwise, the curves in figure 7 would not match. From this analysis, we expect that the resonator finesse is further optimized by minimizing the ROC by decreasing the photoresist d.

Zoom In Zoom Out Reset image size

Our modeling and simulation approach enables us to investigate how the etch times affect the quality of the microcavities, allowing for design optimizations. Consequently, it is important to define objective design criteria through which we evaluate the quality of a microcavity: to that end, firstly, we expect the final opening O to be large enough to contain most of the beam. For this purpose, we require O to be larger than six times $w_\mathrm M$ defined in (7). From a Gaussian beam analysis point of view, this value of O is sufficient to contain three standard deviations of the intensity of the beam inside the resonator. This is enough to obtain resonators with potential finesse larger than the state of the art [10].

Secondly, the best performance of the resonators is obtained when the cavity shape is as parabolic as possible [62]. We model the shape of the microcavity by an even polynomial. We find that a sixth-order polynomial fit $a_2x^2 +a_4x^4 + a_6x^6$, as described in section 4.1, is sufficient, since we seek to manufacture parabolic cavities. Therefore a sixth-order fit suffices to both match the simulated geometry and quantify deviations from the ideal parabolic behavior. We then define our measure of parabolicity P as:

$$\begin{equation} P = \frac{|a_2|}{\sqrt{a_2^2+a_4^2+a_6^2}}. \end{equation} \tag{ 8 }$$

This motivates defining the parabolicity error as $\epsilon_\mathrm{P} = 1 - P$. We set a maximum $\epsilon_\mathrm{P}$ threshold to 10⁻⁶, which is adequate to generate resonators with finesse larger than 10⁶ [60].

Thirdly and lastly, as already discussed in section 4.3, the ROC should be as minimal as possible. Therefore, our optimization targets are summarized as follows:

$$\begin{equation} \left\{\begin{array}{l} O \gt 6w_M, \\ \epsilon_\mathrm{P} \lt 10^{-6}, \textrm{and} \\ \min{(\mathrm{ROC})}. \end{array} \right. \end{equation} \tag{ 9 }$$

The simulation results for the cavity with photoresist $d = 12.4~\mu$m are summarized in figure 8. From figure 8(a)), we see that the cumulative etch time required for the design criterion on the opening O is significantly reduced with increased $t_{\mathrm{resist}}$. This design criterion, shown by a blue arrow in figure 8(a)) will already be fulfilled with the photoresist on at time $t_{\mathrm{resist}} = 600$ s. Combining this insight with the reduction in ROC seen in figure 6, we see good evidence for increasing the $t_{\mathrm{resist}}$ to approximately its double in order to minimize the ROC and achieve a sufficiently large O quickly. However, it must be noted that $t_{\mathrm{resist}}$ cannot be arbitrarily increased due to experimental limitations. Possible issues are photoresist hardening [45, 63] and loss of photoresist structural integrity after significant underetching.

Zoom In Zoom Out Reset image size

Analyzing the $\epsilon_\mathrm{P}$ in figure 8(b)), we highlight the necessity of the second etch step. After the first step, the microcavities have very poor P. It is during the second etch step that we see an improvement of the $\epsilon_\mathrm{P}$. The aforementioned behavior is explained by the second etch step being essentially isotropic without any visibility effects, as illustrated in figure 1(a)). During such an isotropic process, we see that the feature converges to a parabola. This behavior exposes a fundamental trade-off on increasing the second etch time: whereas it improves the $\epsilon_\mathrm P$, it also worsens the ROC. Our analysis suggests that the second etch time should be kept as short as possible, which keeps the $\epsilon_\mathrm{P}$ under a fixed threshold. We thus expect that microcavities of comparable quality are achieved with a reduction of up to $50\%$ of the second etch duration. Combined with the analysis from the first etch step, this results in an overall reduction in etch time of at least $35\%$, representing a significant optimization of microcavity plasma etching processes.