Chapter 6 High Accuracy Hierarchical Grids for Topography Simulation

In Chapter 5, an algorithm is presented that detects areas of a discretized surface with significant geometric variations (i.e., features). As discussed in Chapter 5, typical device topographies originating from process TCAD simulations are composed of extensive areas which are essentially flat and small areas with features. These observations motivate using a simulation domain with varying resolutions (e.g., higher resolutions for features and lower resolutions for non-features) to improve simulation performance.

Locally increasing the resolution of a discretized simulation domain to improve the accuracy of the solution when numerically solving PDEs (e.g., the level-set equation) is a common strategy [19, 20, 21, 22, 23, 24]. Furthermore, these techniques are also employed in other numerical simulations like Monte Carlo simulations [138]. For process TCAD simulations, the foundations are presented in [139]. The technique of locally adapting the resolution of the simulation domain is typically referred to as adaptive mesh refinement (AMR). It should be mentioned here that the term mesh in AMR is a different mathematical object than a conforming mesh (see Definition 2.2.10). To avoid confusion, when the word mesh is used in this work, it refers to Definition 2.2.10, meshes in the context of AMR are called grids.

In level-set based topography simulations that utilize AMR, the entire simulation domain is covered by a so-called base grid (or level 0 grid) and several rectangular refined grids (or level 1 grids) with higher resolutions. Furthermore, it is possible that refined grids (i.e., level 2 grids) with an even finer resolution are nested inside a refined grid (i.e., level 1), which is further discussed in Section 6.1.1. A base grid with all its refined nested grids is called a hierarchical grid. The general term grid refers to all kinds of grids (i.e., the base grid and refined grids). Figure 6.1 illustrates three refinement levels with different resolutions.

A straightforward approach that utilizes hierarchical grids in a level-set based simulation is to cover the entire narrow band around the zero level-set with refined grids [140, 141, 142]. However, this approach still resolves non-features of the geometry with a high resolution, which leads to large areas that do not benefit from the higher resolution. A computationally more efficient way of covering the simulation domain with refined grids is to detect areas of interest in the simulation domain and cover them with refined grids [143]. Therefore, a hierarchical grid placement algorithm is proposed that analyzes the device topography and only covers areas of interest (i.e., features) with refined grids [22, 144].

Figure 6.2 illustrates a level-set function with three features covered by refined grids with higher spatial accuracy (i.e., three times finer).

In order to maximize the performance gains from this hierarchical approach, these refined grids need to be optimally placed, guided by an automatic feature detection method.

The automatic hierarchical grid placement algorithm consists of two parts

• 1. An automatic feature detection step where areas of interest (i.e., features) are marked (flagged) as presented in Chapter 5.

• 2. An automatic grid generation algorithm that clusters the flagged grid points and places refined grids with a finer resolution at these clusters (see Section 6.1.2).

In Section 6.1 sub-grids (i.e., refined grids with additional properties) are defined. Furthermore, nesting criteria, to avoid computational overhead, for these sub-grids and how they are placed in the simulation domain are discussed. Next, in Section 6.2 a hierarchical grid placement algorithm is introduced that places sub-grids based on the detected features of the level-set function (see Section 5.1). This hierarchical grid placement algorithm is then incorporated into the general workflow for topography simulations (see Section 4.5). Section 6.3 evaluates the performance of the hierarchical grid placement algorithm by simulating SEG of silicon-germanium (\(\mathrm {SiGe}\)) fins in narrow oxide trenches.

6.1 Hierarchical Grids

In the introduction to this chapter, the basic terminology for hierarchical grids is introduced (e.g., base grid, refined grids, and level \(n\) grids). For further deliberations in this work, additional terms have to be defined: A positive integer describes the so-called refinement ratio (\(\text {R}_{\text {ratio}}\)) which describes the ratio of the resolution between two grid levels of the simulation domain. When the refinement ratio is even, a grid point on a coarser grid level is covered by the intersection of the \(4\) (2D) or \(8\) (3D) closest grid cells on the next finer level. In the uneven case, the grid point on the coarser level is directly covered by the central grid point of the refined grid with a finer resolution. Figure 6.3 illustrates a grid cell covered by refined grids with different refinement ratios.

The refinement ratio between two grid levels must not necessarily be equal, as suggested in [146]. In the remainder of this thesis, the refinement ratio between grid levels is assumed to be constant with a refinement factor of \(\text {R}_{\text {ratio}}=4\), which is considered an industry standard [147]. However, in Chapter 7, this convention is broken, and the last grid level refinement ratio is dynamically chosen to improve simulation performance.

To solve the level-set equation, at least a star stencil \(\eta _{S}(\mathbf {x})\) is required (see Section 4.1.1). However, points on the border of a refined grid are missing at least one \(\phi \)-value to calculate the finite differences. This, however, is tackled by expanding the definition of a refined grid to a sub-grid [147].

The ghost cells directly neighbor the grid cells at the border of the refined grid (see Figure 6.4).

This construction allows for calculating finite differences even when the finite difference stencil extends beyond the boundary of a refined grid. Furthermore, the ghost cells are used to transfer data between neighboring sub-grids and set boundary conditions. Using ghost cells increases the required memory footprint since each sub-grid has to store its ghost cells. Additionally, if two sub-grids are adjacent, the ghost cells must be stored multiple times. If bigger finite difference stencils are required, the number of ghost cells may be increased. The concept of level \(n\) grids (for \(n > 0\)) can naturally be expanded to sub-grids, e.g., a level 2 sub-grid is a sub-grid with a \(\text {R}_{\text {ratio}}^2\) finer resolution than the base grid.

6.1.1 Nesting

During a typical topography simulation workflow, numerous nested sub-grids are present in the simulation domain. Furthermore, these sub-grids have to exchange data amongst each other (e.g., between sub-grids on the same level or between grid-levels). If the sub-grids could be placed arbitrarily in the simulation domain, a substantial computational overhead would be necessary to handle the data exchange between the sub-grids since all possible relations between sub-grids have to be treated. Thus, the sub-grid placement is restricted to only allow for specific sub-grid configurations. These restrictions allow for efficient handling of the interactions between sub-grids. When a grid \(G_1\) on a grid level overlaps a sub-grid \(G_2\) on a finer level, \(G_1\) is referred to as the parent of \(G_2\), inversely \(G_2\) is referred to as the child of \(G_1\).

The placement restrictions (nesting criteria) for sub-grids used in this work are as follows [147]:

• 1. Sub-grids may not overlap other sub-grids on the same level.

• 2. Each sub-grid has a unique parent grid on the next lower level (i.e., a sub-grid or the base grid).

• 3. Sub-grids have to be aligned to their parent’s grid cells (i.e., each sub-grid is either fully covered by a grid cell of its parent or is not covered at all).

• 4. A sub-grid may not border an area not refined on a lower grid level.

The first two criteria guarantee that each sub-grid has a single unique parent grid, which improves computational efficiency when exchanging data between levels. The last two criteria ensure that there is a gradual refinement from a coarser to a finer resolution.

Figure 6.5 illustrates hierarchical grids that fulfill and do not fulfill the nesting criteria.

6.1.2 Grid Generation

This section discusses how sub-grids are obtained from the flags created by the feature detection algorithm. For the remainder of this section, it is assumed that a single 2D grid with several flagged grid points is given. The extension of the presented algorithm to hierarchical grids is presented in Section 6.2. Furthermore, the general algorithm is independent of the dimension of the grid [148].

The basic grid generation algorithm used in this thesis is a slight variation of the grid generation algorithm of Berger and Rigoutsos [148]. This algorithm considers the flagged and non-flagged grid points in a grid as a black and white image. In this context, the flagged grid points are considered to be black "\(1\)", and the non-flagged grid points are considered to be white "\(0\)". The thus defined image (i.e., the currently examined grid) is called a patch (\(P\)), which is put into a queue. This queue stores all patches until they are considered to be efficient, the efficiency of the patch \(P\) is determined by evaluating

\begin{align} \mathcal {E}(P) = \frac {\text {number of flagged grid points in } P}{\text {total number of grid points in } P}. \end{align} Furthermore, the minimum extent in all Cartesian directions of the patch is determined. If the calculated efficiency value of the patch is smaller than a threshold parameter \(E\) or the minimum extent of the patch has reached a minimum width \(M\), the patch is considered efficient. The efficient patch is removed from the queue and placed into a separate pre-grid list.

If one of the previously discussed criteria is not met, the patch is not efficient and has to be split. The position of the split is determined by calculating the signature (\(\Sigma \)) at the borders of the patch. Starting at each point at the borders of the patch, the number of flagged grid points that lie on a straight line (i.e., a grid line) passing through the patch is calculated. If an element in the signature has a value of \(0\), the patch can be split along this grid line. In case there are no \(0\)s in the signature, an approach related to edge detection in images is used [148]. To that end, the Laplacian (\(\Delta _L\)) of the signature has to be calculated by using the finite difference formula shown in Equation 3.19 with \(\Delta x = 1\). The values of the Laplacian of the signature array are evaluated regarding so-called zero crossings (or inflection points), which are points where the sign of the second derivatives changes. In general, there are multiple zero crossings. Thus, the zero crossing with the biggest numerical change is chosen as the position to split the patch. Furthermore, there may be no clear point at which to split the patch, in this case, the patch is split in half. After the split, all grid lines bordering a new patch with a signature of \(0\) are removed. The newly created patches are added to the queue, and the next patch in the queue is considered.

Figure 6.6 shows an example of how the above discussed grid generation algorithm operates on a 2D grid. In this example, the efficiency parameter is set to \(\mathcal {E}(P) = 0.3\) and the minimum patch side size to \(M = 2\). The signature of the patch has no \(0\)s, so the Laplacian of the signature array has to be calculated. There are two zero crossings (i.e., between -1 and 1 on the x-axis and between -2 and 2 on the y-axis). Thus, the zero crossing on the y-axis is chosen, and the patch is split along the indicated line in Figure 6.6b.

The grid generation algorithm terminates when the queue is empty. Afterwards, for each patch stored in the pre-grid list, a sub-grid, with a finer resolution (according to \(\text {R}_{\text {ref}}\)), is created that covers the patch. The \(\phi \)-values of the newly created sub-grids are calculated in a two-step process. First, the \(\phi \)-values on the newly created sub-grid are initialized using bilinear interpolation from the parent grid; This step initializes the position of the zero level-set on the new sub-grid. Second, a hierarchical re-distancing step is performed that calculates the final \(\phi \)-values of the new sub-grid [149].

6.2 Hierarchical Grid Placement

In this section, the algorithms presented in Section 5.1.2 and Section 6.1.2 are combined into a single algorithm for hierarchical grid placement for level-set functions (see Algorithm 2). Moreover, this algorithm is incorporated into a topography workflow. This algorithm is further referred to as the regridding algorithm. The regridding algorithm starts with a feature detection step that iterates over all grids in the simulation domain and flags the features of the level-set functions (see Section 5.1.2). It should be noted that during the flagging procedure, not only geometric features (i.e., features defined by the curvatures of the geometry) may be of interest. Examples of such non-geometric features are interfaces between multiple material layers (i.e., level-set functions). After the flagging procedure, the automatic grid placement algorithm starts by computing the new sub-grids on the base grid. The flagging and grid generation steps are executed iteratively until a pre-defined number of refinement levels – the general refinement level (\(\text {G}_{\text {ref}}\)) – is reached, which then terminates the grid placement algorithm. If during this process, a sub-grid is found that has not changed position and size from the previous simulation step, the previous data (e.g., the \(\phi \)-values) is copied into the new sub-grid. Furthermore, it is checked if the sub-grids fulfill the nesting criteria (see Section 6.1.1).

The first three criteria are automatically fulfilled due to the construction of the grid generation algorithm (see Section 6.1.2)

• 1. Sub-grids on the same level cannot overlap since each patch is split into smaller non overlapping patches.

• 2. Sub-grids always have a unique parent since the grid generation is executed on each grid on the same grid level individually.

• 3. All sub-grids are aligned to their parents grid cells since only grid cells of the parent grid are refined (see Section 6.1).

Only the fourth criteria has to checked manually. If a sub-grid is found that borders an unrefined area on the previous grid level, the points on the border are flagged as features. Afterwards, the grid generation algorithm is executed again on the previous grid level. The algorithm concludes with a hierarchical re-distancing step [149].

6.2.1 Extended Topography Simulation Workflow

The major difference introduced into the simulation workflow due to the utilization of hierarchical grids is after the advection step, where the regridding algorithm has to be executed. During the advection step, the velocities and the solution to the level-set equation must be calculated on each grid. These steps can be executed on all grids simultaneously. However, the CFL condition (see Equation 4.11) is bound by the grid resolution of the finest sub-grid in the simulation domain. This restriction is necessary since otherwise, the zero level-set on the coarser sub-grids would evolve faster than the zero level-set on the finer sub-grids. Thus, the zero level-set of the level-set function would have two different positions depending on the observed sub-grid. It should be mentioned here that it is possible to evolve the zero level-set on different grid levels with different CFL conditions. However, this has two primary drawbacks. First, the \(\phi \)-values of the coarser grid level have to be interpolated in time to solve the level-set equation on the finer grid level. Second, several edge cases need to be handled, e.g., when another part of the zero level-set evolves into a finer sub-grid during a time step on the coarser level.

Figure 6.7 depicts an updated topography simulation workflow (see Figure 4.7).

The reconstruction of the signed distance function has to be done before the feature detection to guarantee accurate \(\phi \)-values of the level-set function after the advection step, which avoids errors in the feature detection. Furthermore, the regridding does not need to be invoked every simulation time step. Ideally, it should only be invoked when a feature moves out of a sub-grid. Thus, for this work it has been empirically determined, that the regridding procedure should be called every \(4\) time steps.

6.3 Benchmark Example Selective Epitaxial Growth

This section presents an evaluation of the extended topography simulation workflow introduced in Section 6.2.1. By comparing it to simulations of two SEG processes investigated by Jang et al. [137]. To that end, the run-time and accuracy of the simulated device topographies are analyzed. As previously shown in Section 5.2, topographies originating from SEG processes result in geometries with pronounced features which benefit from a higher resolution. Furthermore, in this section, to ensure a well-resolved description of the level-set functions, the material interfaces with the \(\mathrm {SiGe}\) material are considered as features. The simulations were performed with Silvaco’s Victory Process and were executed on the ICS (see Section 4.6.3).

6.3.1 Simulation Setup

In the experiments presented in Jang et al. [137], the authors grow \(\mathrm {SiGe}\) fins inside \(\mathrm {SiO}_{2}\) trenches. The trenches are formed in an initial dry etching step on top of a \(\mathrm {Si}\) substrate. The \(\mathrm {SiGe}\) is grown in a cyclic process inside the \(\mathrm {SiO}_{2}\) trench. \(\mathrm {Si}_2\mathrm {H}_6\) and \(\mathrm {GeH}_4\) is deposited for 15 seconds from flowing source gases. Afterwards, the \(\mathrm {SiO}_{2}\) surface has to be cleaned from nucleated \(\mathrm {Si}_{1-x}\mathrm {Ge}_x\), which is achieved by a 12s etching step by flowing \(\mathrm {Cl}_2\). These two steps compose a single SEG cycle. The composition of the produced \(\mathrm {Si}_{1-x}\mathrm {Ge}_x\) alloy is influenced by the rate of the \(\mathrm {GeH}_4\). This cyclic process leads to the formation of high-quality \(\{1\,0\,0\}\), \(\{1\,1\,1\}\), and \(\{3\,1\,1\}\) crystal facets.

The cyclic SEG process is modeled as a continuous epitaxy process that is perfectly selective (i.e., the deposition rate on the \(\mathrm {SiO}_2\) walls is \(0\)). The process model (i.e., the velocity field) is constructed from experimentally characterized growth rates presented in Jang et al. [137]. Note that the process model does not directly simulate the cyclic process. Nevertheless, the simulated surface of the \(\mathrm {SiGe}\) crystal (that is created after one of these cycles) is referred to as a SEG step or a SEG cycle. One SEG cycle can consist of several level-set simulation time steps (i.e., advection steps).

In Jang et al. [137], the deposition rates on the crystallographic facets for two different alloys \(\mathrm {Si_{0.64}Ge_{0.36}}\) (SEG1) and \(\mathrm {Si_{0.45}Ge_{0.55}}\) (SEG2) were measured. For the \(\{1\,0\,0\}\), \(\{1\,1\,1\}\), and \(\{3\,1\,1\}\) facets (the rates used for the crystal orientation-dependent process model are reported in Table 6.1).

To apply these rates during a topography simulation, the process model utilizes a four-rate Hubbard interpolation to calculate velocity values for all grid points [136].

The epitaxial growth simulations of the \(\mathrm {SiGe}\) crystals SEG1 and SEG2 are performed with three different grid settings, shown in Table 6.2:

A singular grid without sub-grids is utilized when performing the SEG simulations with the Fine and Coarse grid settings. The simulation domains are comprised of \(\approx 3 \cdot 10^5\) grid points when using Fine grid settings and \(\approx 1.9 \cdot 10^4\) grid points when using the Coarse grid settings, which is why it was omitted in the here presented analysis. These two grid settings are intended to create a frame of reference against which the performance and accuracy of the hierarchical approach can be compared. For the sake of completeness, a grid resolution of \(\SI {0.00125}{\micro \meter }\) is considered, which is the grid resolution between Fine and Coarse. However, a simulation using a single grid and a grid resolution of \(\SI {0.00125}{\micro \meter }\) does only perform marginally better than a simulation utilizing the Coarse grid-setting.

The SEG simulations using the Multi-Grid grid settings utilize a base grid with the same spatial resolutions as Coarse grid settings and a general refinement level \(\text {G}_{\text {ref}}\) of \(1\), which leads to the generation of a plethora of sub-grids with the same spatial resolution as the Fine grid settings. On average, a simulation run with Multi-Grid grid settings consists of \(\approx 2.7 \cdot 10^4\) grid points (i.e., the sum of the grid points of the base grid and sub-grids). The interfaces of \(\mathrm {SiO_2}\) and \(\mathrm {SiGe}\) are considered as features (see Section 6.2). For the grid placement algorithm, the efficiency parameter is set to \(\mathcal {E}(P) = 0.7\) and the minimal patch width is set to \(M = 6\).

In what follows, the simulation results attained with the Multi-Grid grid settings are compared to the measurement data reported by Jang et al. [137]. Furthermore, the final \(\mathrm {SiGe}\) crystal surfaces generated with Fine, Coarse, and Multi-Grid grid settings are analyzed and compared against each other. The \(\mathrm {SiGe}\) crystal surface generated with the Fine grid settings is considered the comparison’s reference surface since it has the highest numerical accuracy. Therefore, the error introduced by the Coarse and Multi-Grid grid settings is investigated by calculating the smallest Euclidean distance (i.e., \(d_E\)) between each surface point of the \(\mathrm {SiGe}\) crystal simulated with Coarse and Multi-Grid grid settings and the closest point on the reference surface (\(d_E\)-error).

6.3.2 Example 1

The simulated \(\mathrm {SiGe}\) surfaces after several SEG cycles using the SEG1 process model (see Table 6.1) and the Multi-Grid grid settings compared to the experiment are shown in Figure 6.8.

The simulation results are in good agreement with the experiment after all three measured cycles. Figure 6.9 shows the final surfaces (i.e., after 47 SEG cycles) of the SEG1 simulation with all grid settings reported in Table 6.2.

The Fine and Multi-Grid simulation results are in good agreement. On the other hand, the simulated \(\mathrm {SiGe}\) crystal surface utilizing the Coarse grid settings does not match the other simulated surfaces. The error in the crystal surface is most pronounced in the mismatch of the peak of the \(\mathrm {SiGe}\) crystal.

The \(d_E\)-error between the \(\mathrm {SiGe}\) crystal surfaces is reported in Figure 6.10. As again indicated by the spike in the \(d_E\)-error, the simulation results obtained with the Coarse grid settings do not offer sufficient accuracy compared to the other grid settings. Furthermore, the \(d_E\)-error between the Multi-Grid and Fine grid settings is negligible (i.e., it is much smaller than the grid resolution), showing the anticipated increase in accuracy of the hierarchical approach. The \(d_E\)-error of the surface using the Multi-Grid grid settings is as big as the expected error from a simulation run using the Fine grid settings.

The simulation run-times of the SEG1 and SEG2 simulations run with the three studied grid settings are shown in Table 6.3:

The results show that the Fine grid settings have the worst simulation run-times. This is expected since the high resolution used in the entire simulation domain unnecessarily resolves non-features that do not benefit from the high resolution. In contrast and as expected, the simulation using the Coarse grid-setting results in the fastest run-times. However, as discussed previously, the quality of the final simulated surface is insufficient. The Multi-Grid grid settings improves the simulation performance by 32% compared to the Fine grid settings.

The flagged grid points and thus generated sub-grids after 24 SEG cycles, using the SEG1 process model, are shown in Figure 6.11.

During this simulation step of the SEG process, a new feature forms between the \(\{1\,1\,1\}\) and \(\{3\,1\,1\}\) crystal facets, which are detected by the feature detection algorithm. The hierarchical grid placement algorithm places the sub-grids accordingly, as indicated by (a) in Figure 6.11.

6.3.3 Example 2

The simulated \(\mathrm {SiGe}\) surfaces utilizing the SEG2 process model and Multi-Grid grid settings compared to the experiment are shown in Figure 6.12.

The final surfaces (i.e., after 55 SEG cycles) of the \(\mathrm {SiGe}\) crystal using all grid settings given in Table 6.2 are shown in Figure 6.13.

The error between the Coarse and Fine grid settings is not as pronounced as in the first experiment (i.e., SEG1), although there is still a mismatch between the surfaces, especially at the peak of the crystal. The surfaces obtained from the Multi-Grid and Fine simulations are again in good agreement.

Figure 6.14 reports the \(d_E\)-error of the crystal surfaces simulated with the SEG2 process model.

The \(d_E\)-error of the Coarse grid settings is again larger than the error of the Multi-Grid grid settings. Furthermore, the minimal \(d_E\)-error between the surfaces is bigger than in the SEG1 experiment. These results again confirm the anticipated increase in accuracy obtained by the hierarchical grid with a negligible \(d_E\)-error.

Considering the run-times reported in Table 6.3 for the SEG2 process model, the observations made in the previous experiment still hold. The Coarse grid-setting again has the fastest run-time, with the drawback of an inaccurate description of the \(\mathrm {SiGe}\) crystal surface. In the SEG2 experiment the Multi-Grid grid settings improve the simulation run-time by 58% compared to the Fine grid settings.

Another point that requires discussion is the difference in total run-time between the SEG1 and SEG2 experiments (see Table 6.3). The SEG2 experiment simulates 55 SEG cycles, while the SEG1 experiment only simulates 47 cycles with more than double the simulation run-time compared to SEG2. This difference in run-time is explained by the complexity of the velocity field (i.e., process model), which influences how far the zero level-set can propagate in a single time step. A closer investigation of the actual time steps the level-set method has to calculate during the SEG simulation shows that for the Multi-Grid grid settings the SEG1 simulation has to perform \(2259\) time steps compared to \(819\) for SEG2. This explains the difference in simulation run-times between SEG1 and SEG2, which is due to the different number of simulation time steps that must be performed.

Figure 6.15 shows the flagged grid points and sub-grids generated by the hierarchical grid placement algorithm after 33 SEG cycles using the SEG2 process model.

The process model for the SEG2 process does not generate a new feature during the SEG process. Thus no additional sub-grids are needed on the crystalline surface. Still, the peak of the \(\mathrm {SiGe}\) crystal is detected as a feature, and a sub-grid is correctly placed accordingly.

6.4 Summary

The feature detection algorithm presented in Chapter 5 has been combined with an AMR algorithm to create a hierarchical grid placement algorithm. It has been discussed how the AMR algorithm determines the size and position of sub-grids. Furthermore, nesting criteria for nested sub-grids as well as a strategy to calculate finite differences at the borders of these sub-grids, have been discussed. The hierarchical grid placement algorithm has been integrated into the topography simulation workflow. To demonstrate the efficiency of the hierarchical grid placement algorithm for topography simulations, representative and practically relevant simulations of selectively grown epitaxial \(\mathrm {SiGe}\) fins in oxide trenches have been performed. The SEG process with two different parameter sets has been successfully simulated. During the simulation of the SEG processes, the hierarchical grid placement algorithm allows for a low base grid resolution combined with sub-grids with a higher resolution which are optimally placed to cover the features of the \(\mathrm {SiGe}\) crystal. The hierarchical grid placement algorithm presented in this chapter improves the simulation run-time by 32% and 58% with respect to the two process models, SEG1 and SEG2, respectively compared to a simulation run with a single high resolution grid while also upholding accurately.