Chapter 4 The Level-Set Method

This chapter starts with the theoretical (mathematical) background of the level-set method (cf. Section 4.1), clearly defining the used terms such as level-set function and signed-distance function. In Section 4.2, the reference implementation of the level-set method within the context of a simulation tool is presented, providing a detailed discussion of all the computational steps in dedicated subsections. The computational steps range from the creation of the level-set function representing the interface and the movement or deformation of said interfaces, and finally to the extraction of an explicit representation of the interfaces from the full simulation domain, which is then used for further process or device TCAD simulations. The chapter concludes with a short overview of the available process TCAD simulators which use the level-set method (cf. Section 4.3).

4.1 Theoretical Background

This section discusses the implicit representation of geometries using a function and their movement (deformation).

4.1.1 Level-Set Function

For a given \(\Omega \subset \mathbb {R}^n\) (process TCAD simulations typically consider \(n=2\) or \(n=3\)) and an interface \(\Gamma \subset \Omega \) separating the inside \(\Omega ^-\subset \Omega \) from the outside \(\Omega ^+=\Omega \setminus \Omega ^-\). The level-set method describes the interface separating the inside from the outside implicitly by a function \(\Phi \left (\vec {x}\right )\). This function also known as the level-set function, has to be continuous and fulfill

\begin{align} \Phi \left (\vec {x}\right )\begin{cases} <0 & \quad \text {for } \vec {x}\in \Omega ^-, \\ =0 & \quad \text {for } \vec {x} \in \Gamma , \\ >0 & \quad \text {for } \vec {x}\in \Omega ^+. \end {cases}\label {eq::lsf} \end{align} Because this work only considers the 2D and 3D case, let \(\vec {x}=\left (x,y,z\right )\), where the z-component is only considered in the 3D case. Figure 4.1 shows a 2D example of an interface. The choice on which side of the interface the sign is negative is arbitrary. The convention in this thesis is that the inside shall have the negative sign. Also, the choice to select the zero-level-set as the interface is arbitrary (any other value is possible as well), but by choosing zero the sign of the function is sufficient to discriminate between inside and outside.

The gradient of \(\Phi \), if it exists (on corners and edges of \(\Gamma \) the gradient does not exist), is given by

\begin{align} \nabla \Phi =\left (\frac {\partial \Phi }{\partial x},\frac {\partial \Phi }{\partial y},\frac {\partial \Phi }{\partial z}\right ), \end{align} and is perpendicular to the iso-contours of \(\Phi \), including the special iso-contour for the value zero (zero-level-set) of \(\Phi \). Therefore, normalizing the gradient yields the normal vector¹ of the interface

\begin{align} \vec {n}=\frac {\nabla \Phi }{|\nabla \Phi |}. \end{align} The so defined normal vector points outwards because the gradient points in the direction of increasing \(\Phi \). \(\Phi \) is by definition (4.1) smaller on the inside than on the outside. The definition of the normal vector via the normalized gradient allows for a straightforward embedding of the interface normal vector \(\vec {n}\) in \(\Omega \). Therefore opening possibilities for a geometric interpretation of the level-set function off of the interface, e.g., curvature in the entire domain \(\Omega \). Unfortunately, the minimal requirements to the level-set function do not allow for a geometrical interpretation of the normal vector at points which are not on the interface \(\Gamma \). Thus, practical level-set simulations typically require more properties of the level-set function, discussed in the next section.

4.1.2 Signed-Distance Function

In practical applications the level-set function \(\Phi \) is often chosen as a signed-distance function relative to \(\Gamma \)

\begin{align} \Phi \left (\vec {x}\right )=\begin{cases} -d(\vec {x},\Gamma ) & \quad \text {for } \vec {x}\in \Omega ^-, \\ \hfill 0 & \quad \text {for } \vec {x} \in \Gamma , \\ +d(\vec {x},\Gamma ) & \quad \text {for } \vec {x}\in \Omega ^+, \end {cases}\label {eq::sdf} \end{align} with \(d\) the Euclidean distance to the interface \(\Gamma \). The Euclidean distance between two points is given by

\begin{align} d\left (\vec {x},\vec {y}\right )=\sqrt {\sum _{i\in \{x,y,z\}}(\vec {x}_i-\vec {y}_i)^2}, \end{align} with \(\vec {x}_i\) the component of the vector \(\vec {x}\) in the corresponding spatial direction. The distance between a point and a set (i.e., the interface \(\Gamma \)) is given by the infimum of the Euclidean distance over all points of the set

\begin{align} d\left (\vec {x},\Gamma \right )=\inf _{\vec {y}\in \Gamma }d(\vec {x},\vec {y}). \end{align}

The choice, to use a signed-distance function, allows for geometrical interpretation of the normal even for points not on the interface. The normal points away from the closest point on the interface on the outside, and to the closest point on the interface on the inside. If \(\Phi \) is a singed-distance function \(|\nabla \Phi |=1\) holds.

This fact is easily reasoned by considering a point \(\vec {x}\) and the corresponding closest point on the interface \(\vec {x_\Gamma }\). All points \(\vec {y}\) on the shortest path connecting \(\vec {x}\) and \(\vec {x_\Gamma }\) have the same point \(\vec {x_\Gamma }\) as the closest point on the interface. This is due to the triangle inequality. Thus the path connecting \(\vec {x}\) and \(\vec {x_\Gamma }\) is the path of steepest descent for the level-set function \(\Phi \), evaluating to \(-\nabla \Phi \). Furthermore, because \(\Phi \) is scaled according to the Euclidean distance it follows that \(|\nabla \Phi |=1\).

The representation using a signed-distance function allows for a faster computation of the normal vector of the interface, because the normalization step is dispensable

\begin{align} \vec {n}=\frac {\nabla \Phi }{\underbrace {|\nabla \Phi |}_{=1}}=\nabla \Phi . \end{align}

In Figure 4.2 two different level-set functions for the same interface (a circle) are shown. Figure 4.2a shows an arbitrary level-set function whereas Figure 4.2b shows the use of a signed-distance function. Figure 4.2c and Figure 4.2d show the corresponding iso-contours. In case of the signed-distance function the iso-contours are evenly spaced, whereas in the other case they tend to cramp up or spread out. If the level-set method is solved analytically, this is not an issue, but in case of a discretization and numerical procedures issues arise. Numerical issues arising from such steep or flat gradients of the level-set function are:

• • Steep gradients may exceed the numerical representation of floating-point numbers.

• • Flat gradients are prone to distortions of the interface position (small perturbations of the level-set function lead to enormous perturbations of the interface position).

One point in the domain for the signed-distance function in Figure 4.2 is special, i.e., the center of the circle (the apex of the cone), because it is the only point for which \(|\nabla \Phi |\) is undefined. Generally, for a signed-distance function the gradient for all points on the skeleton of \(\Gamma \) is undefined. The skeleton of an interface \(\Gamma \) consists of all points \(\vec {p}\in \Omega \) which have more than one closest point on \(\Gamma \), e.g., for a circle the center or for a square the diagonals [106].

This seems to be a critical drawback when using a signed-distance function, but as the equations under consideration in this work are generally true, e.g., \(|\nabla \Phi |=1\) holds almost everywhere (except for a negligible subset of \(\Omega \)).

The advantages of the geometric interpretation of the level-set function and its numerical robust gradient outweigh the drawback not being able to define the gradient everywhere. An almost everywhere true equation, e.g., \(|\nabla \Phi |=1\), may still be approximated numerically, if the approximation ’Fails in a graceful way’ [18], meaning that the failure does not cause a deterioration of the underlying numerical method.

4.1.3 Interface Movement

Assume the velocity \(\vec {V}\) is given for each point of the interface \(\Gamma \). The movement of the interface is given by the movement of all interface points, whereas the movement for a single interface point \(\vec {x}\) is described by

\begin{align} \frac {d\vec {x}}{dt}=\vec {V}\left (\vec {x},t\right ).\label {eq:lagrange} \end{align} Such a description of the interface movement would require an explicit representation of the interface, because the movement of individual points of the interface is described. The level-set method is an implicit approach, therefore, the approach has to be adapted. The level-set function \(\Phi \) is used in the description of the interface movement. For the interface movement a time dependency is introduced to the level-set function \(\Phi \left (\vec {x},t\right )\). The movement of the interface is described by the advection (convection) equation

\begin{equation} \frac {\partial \Phi \left (\vec {x},t\right )}{\partial t}=\vec {V}\left (\vec {x},t\right )\cdot \nabla \Phi \left (\vec {x},t\right ),\label {eq:ls} \end{equation}

where \(\Phi \) is the well-known signed-distance function and \(\vec {V}\) is a velocity field describing the interface movement. Equation (4.9) is also known as the level-set equation. The interface is moved because the zero-level-set of \(\Phi \) changes over time.

In contrast to the formulation in (4.8), which presents the movement from a Lagrangian specification (the observer follows points), (4.9) uses the Eulerian specification (observer watches which points pass by). Both specifications are related via the equation

\begin{align} \Phi \left ( \vec {X} \left (\vec {x}_0,t \right ),t\right )=\frac {\partial \vec {X}}{\partial t}\left (\vec {x}_0,t \right ), \end{align} where \(\vec {X}\left (\vec {x},t\right )\) is the position of the point \(\vec {x}\) according to (4.8) at time \(t\) and \(\vec {x}_0\) is a generic point on the interface.

In cases where strictly only the interface is of interest a scalar velocity field is sufficient to prescribe the interface movement. The interface position is only affected by the to the interface orthogonal component of the velocity field. Let \(\vec {V}=v\vec {n} + \vec {P}\), where \(\nabla \Phi \perp \vec {P}\) holds and \(v\) is the scalar velocity in normal direction of the interface. Then (4.9) simplifies to

\begin{align} \frac {\partial \Phi }{\partial t} & =\vec {V}\cdot \nabla \Phi , \\ & =\left (v\vec {n} + \vec {P}\right ) \cdot \nabla \Phi , \\ & =v\underbrace {\vec {n} \cdot \nabla \Phi }_{=|\nabla \Phi |} + \underbrace {\vec {P} \cdot \nabla \Phi }_{=0}. \end{align} Additionally, in case of a signed-distance function, (4.9) simplifies further to

\begin{equation} \frac {\partial \Phi }{\partial t}= v.\label {eq:ls1} \end{equation}

This final formulation allows for a short and elegant way to describe the interface movement for a given velocity field.

In the next sections the computational steps of the used reference simulator are presented.

4.2 Reference Simulation Workflow

This section provides details of the considered reference process TCAD simulation workflow which is based on the level-set method. In essence, the simulation workflow consists of three main parts (cf. Figure 4.3):

• • The initialization, in which the level-set functions are set up (Section 4.2.1).

• • The main time loop (Section 4.2.2 – Section 4.2.7) where the device structure (topography) is advanced in small time steps (solving the level-set equation (4.9) in time steps of size \(\Delta _t\)).

• • The finalization in which an explicit representation of the material regions is extracted (Section 4.2.8).

The next sections provide a more in-depth overview of the individual steps shown in Figure 4.3.

4.2.1 Initial Interfaces

The first step of a level-set simulation is to create the level-set functions to represent the material regions. The spatial discretization approach considered here (hierarchical grids) allows to discuss the implementation as if it would use as single Cartesian grid, except for Re-Gridding which adapts the blocks on the hierarchical grid in position and size to the changing material regions.

Figure 4.4 shows a discretized level-set function on a Cartesian grid. On each grid point an approximation of the signed-distance to the interface is stored. Each grid point is located in the center of a cell (black square). Such a level-set function may be used to describe a single material region. Technically, a level-set function describes two regions: 1) A material region and 2) the complement to the material region, e.g., the device structure and the void above the device structure in the simulation domain.

The discuss continues with a detailed description of material representation approaches.

Material Representation

As shown by the thermal oxidation example (cf. Section 1.1) device topographies consist of more than a single material region. A level-set function may only model a single material interface, thus several level-set functions are necessary to represent multiple material regions (cf. Figure 4.5).

There are different approaches how the level-set functions are configured to achieve this goal. The straightforward approach using a dedicated level-set function encapsulating the material region for each material as shown in Figure 4.5a, has some drawbacks which are further illustrated in Figure 4.6. At the interface between two material regions non-physical voids may form (due to numerical inaccuracies of the level-set functions), because two level-set functions represent the same material interface (cf. Figure 4.6a). Sometimes no voids materialize, but the material regions defined by the level-set functions overlap. However, typical process TCAD simulations do not consider (allow) alloys, i.e., the material has to be unique at each point (cf. Figure 4.6b).

Also triple-junctions (points where three material regions meet) are destined to form non-physical voids, because a numerical implementation causes slight unavoidable rounding of corners (cf. Figure 4.6c). The rounding is related to the grid resolution, if the grid resolution approaches zero the rounding vanishes completely.

Several strategies to avoid the aforementioned voids were developed, primarily driven by research on multiphase flows [107, 108, 109, 110]. The material-specific level-set functions are either held together, e.g., deliberately changing the velocity field before the Advection to avoid the forming of voids or forced together, e.g., using Boolean operations after the Advection to remove overlap of the level-set functions.

In general, the focus of the various research strategies is on triple-junctions in 2D fluid dynamic simulations. The previously discussed approaches all require \(M-1\) level-set functions for \(M\) different materials. A level-set function separates two material regions, thus regions belonging to no level-set function are not explicitly represented. In context of process TCAD simulations the void (vacuum or gas) above the to-be-simulated structure is typically the material that is not explicitly represented.

In [111] a conceptually different approach is presented, which uses for each material pair interface a dedicated level-set function leading to a maximum number of level-set functions of \({M(M-1)}/{2}\) for \(M\) materials and all materials having a pairwise interface. A sophisticated voting mechanism decides which of all these level-set functions is used for the actual computations. This enables to represent triple-junctions without voids and overlaps. However, the high number of level-set functions required in this approach represents a computational burden. The most widely-used approach in process TCAD simulations is to use additive level-sets [3, 20] (cf. Figure 4.5b). Instead of storing individual level-set functions representing a dedicated material, the level-set functions store a union of materials. The approach is considered in this work and further presented in the following paragraphs.

Material Representation in Process TCAD Simulations

In an additive level-set approach the first material (used in a simulation) is represented with a single level-set function, as described in Section 4.1.1. Every time a new material is added (can happen several times during a practical simulation workflow) a new level-set function describing the union of all previously present material regions and the new material is added. The defined level-set function always represents the interface between the structure (device) and the vacuum or gas region.

The additive level-set approach prevents the formation of non-physical voids, because at a material boundary only a single level-set function defines the interface, in contrast to the straightforward approach (cf. Figure 4.5a), where two level-sets are present (one for each material).

Additionally, the additive level-set approach allows representing material regions with thicknesses less than the grid resolution. The union with the other materials creates a thicker material region which can be represented by a level-set function. The straightforward approach representing material regions individually is not able to reliably store such thin material regions: Storing such thin material regions is only possible, if the material region is aligned to the grid and has a symmetric offset to the grid points it encloses. Thus, only planar structures are possible in the straightforward approach, if thin layers are considered.

The additive level-set approach allows simulating physically accurate etching processes where thin etch stop layers (a thin material region which is hardly affected by the etching process) are employed, without the necessity for unfeasible high spatial discretization. The material regions of the etch stop layers may then be reconstructed via Boolean operations after an explicit interface representation has been extracted (explicit interfaces are not bound to a grid and therefore have no restriction originating from the grid resolution on their thickness).

A modification to the additive level-set approach to simulate deposition processes, where the process model is highly dependent on the interface normals (e.g., epitaxial crystal growth), is derived in [112]. There, a different strategy to unionize the material regions is used, with the goal that the outer most level-set, which is in this case not the wafer surface, has the correct local curvature (convexity and concavity) enabling the formation of crystal facets. However, Boolean operations allow a straightforward conversion to the additive level-set approach.

The additive level-set approach heavily relies on Boolean operations, which are straightforward for level-set functions and discussed in the following.

Boolean Operations

Boolean operations on level-set functions allow for efficient and high performing implementations [106], because their computation is based on a point-wise evaluation of minimum, maximum or multiplication by \(-1\). For two level-set functions \(\Phi _1\) and \(\Phi _2\) the operations are defined by

\begin{align} \Phi _1 \cup \Phi _2 & = \min \left (\Phi _1, \Phi _2\right ),\label {eq:union} \\ \Phi _1 \cap \Phi _2 & = \max \left (\Phi _1, \Phi _2\right ),\label {eq:intersection} \\ \Phi _1 \setminus \Phi _2 & = \max \left (\Phi _1, -\Phi _2\right ),\label {eq:difference} \\ \Phi _1^c & = -\Phi _1.\label {eq:complement} \end{align} The union of those two level-set functions is given by (4.15), the intersection by (4.16), the difference by (4.17), and the complement of a single level-set function by (4.18). Boolean operations do not preserve the signed-distance property of the level-set function.

4.2.2 Process Model

The process model is fully responsible to define the changes to the topography by determining how the interfaces shall be transformed. In other words, the process model captures all the physics behind the simulation.

The complexity of the chosen process model varies heavily depending on the required accuracy and simulated process. The least complex process models are Boolean operations (see Section 4.2.1) and uniform deposition or etching models (equivalent to the computation of a constant offset of the zero-level-set). More complex process models are interface normal dependent models [113, 114, 115, 116] modeling anisotropic etching processes or epitaxial growth processes. Material flow processes, such as thermal oxidation [15], require the solution of a Navier-Stokes equation. There are also process models which use visibility calculations (ray tracing) to simulate direct particle transport [117, 118, 20, 13]. More sophisticated models extend their particle transport process models further to also account for an additional external flow [26, 119, 120]. The particle transport is essential for so-called reactive ion etch processes, in which the wafer surface is bombarded with ions, kinetically removing material and consequently realizing the creation of high aspect ratio devices.

The usage of the level-set method strictly decouples the process model from the interface advection. This decoupling enables a straightforward switching of the process model, while still using the same material representation and interface evolution. Therefore, comparing different process models of varying physical complexity and accuracy for the same process step is viable, which is important for practical process TCAD simulations to fine-tune predictions of fabrication processes.

In the scope of the reference simulation workflow the interaction of the process model and the level-set method uses a standardized API, which is specified in the next section.

4.2.3 Interface Velocity

The API between the process model and the level-set method is realized via Cross Points. Cross Points are intersections of the interface with grid lines (cf. Figure 4.7). Cross Points are chosen because they are the minimal requirement to a process model: A process model must be able to describe the movement of points on the interface. Some process models (e.g., isotropic deposition) directly yield velocities for all grid points, whilst others (e.g., models using visibility calculation) are not able to provide meaningful velocities at grid points off the interface.

A velocity given only at Cross Points would enable solely a Lagrangian movement of the interface (cf. Section 4.1.3). Thus, the velocity of the interface has to be extended to enable the Eulerian movement necessary for the implicit interface representation employed by the level-set method. This extension step is called Velocity Extension and is discussed in the following.

4.2.4 Velocity Extension

The computational step Velocity Extension is used to extend the velocity from the Cross Points to the grid points of the computational domain. This is necessary because the solution of (4.9) requires the velocity to be defined not only directly at the interface but also on the entire computational domain. The requirement for this extended velocity is to describe the interface movement, therefore, at the interface the extended velocity field and the interface velocity should match. If the extended velocity and the interface velocity do not match at the interface, they describe different interface movements.

To that end, Velocity Extension assigns each grid point the velocity of the closest interface point [121]. A direct computation of the closest interface point and the corresponding Cross Point on an implicit interface is prone to numerical errors. The closest point on the interface \(\vec {x}_{\Gamma }\) of \(\vec {x}\) is computed by \(\vec {x}_{\Gamma }=\vec {x}-\Phi (\vec {x})\nabla \Phi (\vec {x})\) (this holds only for a signed-distance function). Even minor deviations of \(\Phi \) from a signed-distance function are able to break the approach. Especially, near corners of the interface (on different sides of a corner the velocities may differ by a large margin, e.g., caused by a process model using visibility calculations and one side of the corner is shadowed) and for points further away form the interface this approach becomes very inaccurate due to numerical issues (the uncertainty of the computed point scales with the distance to the interface).

An alternative to direct computation of the closest interface point is an extension from the interface constant along the normal direction of the interface, in an outwards marching manner. The approach is formulated by a PDE (a boundary value problem) known as velocity extension equation [122]

\begin{align} \nabla \Phi \left (\vec {x},t\right ) \cdot \nabla V \left (\vec {x},t\right ) & =0,\quad \vec {x}\in \Omega ,\label {eq:velext} \\ V \left (\vec {x},t\right ) & =V_I \left (\vec {x},t\right ),\quad \vec {x}\in \Gamma . \end{align} The given velocity on the Cross Points (interface) is denoted by \(V_I\) whilst the extended velocity, which is used for the advection is denoted by \(V\).

The solution to this PDE is constant along the normal direction (\(\nabla \Phi \) approximates the normal direction) of the interface. Because the change of the velocity (expressed by \(\nabla V\)) is orthogonal to the normal direction (4.19) each point gets the velocity of the closest point of the interface assigned by solving above equations. Extending according to (4.19) has additional numerical benefits, such as a reduced distortion of the signed-distance property, which is particularly relevant for Advection, as discussed in the following. The details and proposed algorithms for Velocity Extension (solving the PDE (4.19)) are presented in Chapter 5.

4.2.5 Advection

After Velocity Extension the level-set equation,

\begin{align} \frac {\partial \Phi \left (\vec {x},t\right )}{\partial t}=\underbrace {\vec {V}\cdot \nabla \Phi \left (\vec {x},t\right )}_{H(\vec {x},\Phi ,\nabla \Phi ,t)}, \tag {\ref {eq:ls} revisited} \end{align} is well-defined in the computational domain. The right-hand side \(H\) is called Hamiltonian in the context of Hamilton-Jacobi equations [73, 18]. The level-set equation is discretized in time². For a time \(t^n\), let \(\Phi ^n=\Phi \left (t^n\right )\) be the description of the interface at \(t^n\) and let \(\Phi ^{n+1}=\Phi \left (t^{n+1}\right )\) be the description of the interface after a small time step \(\Delta _t=t^{n+1}-t^n\). A first-order accurate solution of (4.9) is given by the forward Euler method

\begin{align} \frac {\Phi ^{n+1}-\Phi ^n}{\Delta _t} & =H(\vec {x},\Phi ^n,\nabla \Phi ^n,t^n), \end{align} with \(H\) the Hamilton evaluated at time \(t^n\).

Higher order schemes in time are the total variation diminishing (TVD) Runge-Kutta (RK) schemes, introduced in [123] and further developed in [124]. The schemes combine several sequential forward Euler steps effectively canceling low order error terms, thus yielding a higher order in time. Those TVD RK schemes avoid spurious oscillations as long as the underlying forward Euler scheme does not introduce them. However, the numerical benefit of schemes of order four or higher does not contribute significantly to the accuracy of practical simulation problems [18].

For the discretization of the Hamiltonian several schemes have been developed: Lax-Friedrich [125], Godunov [126], and Roe-Fix with entropy correction [123]. They all have in common that they make use of a numerical Hamiltonian \(H'\) which includes artificial dissipation to damp spurious oscillations in the solution. This allows the basic forward Euler step to be stable.

The typically used Lax-Friedrich scheme computes the artificial dissipation globally giving high dissipation, but often lead to smoothed solutions. High dissipation is not desirable, because the level-set method should not affect the process model, e.g., smoothed structures. To counter the high dissipation, schemes which compute the dissipation locally, like the Local-Lax-Friedrich scheme [127], were developed. A recently developed scheme for anisotropic etching process TCAD simulation considers more grid points than the the Local-Lax-Friedrich scheme to determine the artificial dissipation [115, 112].

The numerical solution of the forward Euler method (4.9) is explicit, thus the limit on the size of the time step is only given by the Courant-Friedrichs-Lewy (CFL) condition

\begin{align} \Delta _t<\frac {\Delta _{\vec {x}}}{\max {\left |\vec {V}\right |}}, \end{align} where the maximum over the computational domain is chosen and with \(\Delta _{\vec {x}}\) a measure for the grid resolution, e.g., the maximum of the grid resolution in all spatial dimensions. Often the signed-distance property of the level-set function is lost during the advection [128]. This is often the result of a bad velocity function (i.e., a velocity not fulfilling (4.19)) or the result of accumulated numerical errors.

4.2.6 Re-Distancing

Re-Distancing restores the signed-distance property of a level-set function. This is necessary because the signed-distance property is distorted by the previous computational step, the Advection. Re-Distancing avoids numerical issues arising from steep and flat gradients of \(\Phi \) (cf. Section 4.1.2).

The approaches considered here to compute the signed-distance stem from an generalized mathematical problem which is known as the Eikonal equation [129]

\begin{align} |\nabla \Phi \left (\vec {x}\right )| & =F\left (\vec {x}\right ) \quad \vec {x} \in \Omega ,\label {eq:eikonal} \\ \Phi \left (\vec {x}\right ) & =G\left (\vec {x}\right ) \quad \vec {x}\in \Gamma .\label {eq:eikonal_2} \end{align} The Eikonal equation describes a wave front emerging from \(\Gamma \) and marching through \(\Omega \). The wave speed is given by \(\frac {1}{F\left (\vec {x}\right )}\), which has to be positive to be well-defined (negative or zero wave speed would not allow the wave front to reach the entire \(\Omega \)).

The solution to \(\Phi \) in this context describes the shortest travel time for the wave emerging from \(\Gamma \). The initial values given by \(G\) on \(\Gamma \) determine the departure time from \(\Gamma \). Early departures are given by \(G<0\) and late departures by \(G>0\).

If the wave speed is uniformly equal to one and the departure time is zero, the travel time is equal to the traveled distance. Thus, for Re-Distancing the PDE

\begin{align} |\nabla \Phi \left (\vec {x}\right )| & =1 \quad \vec {x} \in \Omega ,\label {eq:redist} \\ \Phi \left (\vec {x}\right ) & =0 \quad \vec {x}\in \Gamma , \end{align} is solved. Note, this would just compute the distance field (not signed-distance as both sides of the interface are positive), thus the computed values on \(\Omega ^-\) have to be multiplied by -1 to get the signed-distance function. The algorithmic details, the implementation details, and the developed advanced parallelization strategies are presented in Chapter 6.

4.2.7 Re-Gridding

Re-Gridding adapts the blocks on the hierarchical grid to fit the new interface position (regions requiring fine spatial resolution). Re-Gridding is split into two sub-steps:

• • Flagging of the regions which need a higher spatial resolution.

• • Clustering those regions (covering with blocks) in order to create the hierarchical grid.

Figure 4.8 shows the interface displacement by the advection and the two sub-steps of Re-Gridding.

Flagging

The flagging sub-step selects grid points based on the distance to the interface (approximated by their level-set value \(\Phi \)), the difference of normals on grid points to normals on neighboring grid points (detecting regions with curvature, i.e., corners and edges), and distance to other level-set functions present in the simulation domain (detection of material borders and triple-junctions). Flagging of regions with high curvature is particularly important, because the maximum representable curvature is indirectly proportional to the spatial resolution [18]. This is essential because the often present sharp corners in process TCAD simulations would otherwise become artificially rounded.

Clustering

Clustering is important to efficiently transform the flagged points to actual blocks of the hierarchical grid. The main challenges are to create as few blocks as possible containing the least amount of grid points, whilst still covering all flagged grid points. The run-time of a step in the main time loop almost linearly depends on the number of discretized points, thus minimizing them is important for the overall run-time efficiency of the overall simulation.

This concludes the computational steps of the main time loop, which are repeated until the desired end time of the simulation is reached. After all time steps are done the simulation results are ready to be transferred to the next process step of the manufacturing of the device. This could be starting a new time loop with a different process model or in case the device manufacturing is completed the topography is extracted for the subsequent device TCAD simulation. This latter step requires the interface to be extracted, as is discussed in the following.

4.2.8 Interface Extraction

Device TCAD simulations require an explicit description of the materials regions, thus the implicit interfaces are extracted to an explicit description of the material regions. The preferred explicit representation of the material regions for device TCAD simulations is a polygonal surface representation (e.g., triangle mesh). The polygonal representation is reached by computing an explicit interface from the implicit representation of the level-set function and then use Boolean operations to restore the material regions from the additive level-sets. The dominant algorithm for explicit interface extraction is the marching cubes algorithm [130] in three dimensions, the analog algorithm in two dimensions is called marching squares. Other approaches are for example the cut cell method [131], which is conceptually similar to the marching cubes algorithm, but is tailored towards fluid dynamic simulations. Because this thesis is not focused on this computational step the reader is referred to [132] for an comprehensive overview of advancements to the base marching cube algorithm presented in the following.

Starting from a Cartesian grid, chunks of eight neighbor points (corners of a cube) are used to determine the polygonal representation of the interface passing through the cube. The sign of the points determines which of the \(2^8\) possible polygon templates is chosen. In Figure 4.9 eight different sign combinations with their corresponding polygon templates are shown. The number of templates is reduced through symmetry exploitations (reflection, rotation, and mirror) to a minimum of 14 cases. If not all symmetries are exploited the number of cases is higher. The explicit vertices of the polygons are given by the zero-crossings along the edges of the cube (they are identical to the Cross Points introduced earlier). In a last step, the polygons of all chunks are merged together forming the polygonal representation of the full interface. The so-created polygonal surface representations are often of poor quality (e.g., triangles with high aspect ratios, neighboring triangles with significant different diameters).

Therefore, the initially marching-cubes generated explicit polygonal surface representations are typically optimized and consequently volume-meshed to enable the full-range of subsequent device TCAD simulations (e.g., simulations of charge carrier densities under bias conditions inside the device and at device contacts) [133].

In the next and final section of this chapter, an overview of software tools is given, which implement level-set based process TCAD simulations.

4.3 Software

The advantages of the level-set method in tracking topography evolution lead to the development of several commercial and open source software tools, referred to as process TCAD simulators. They utilize the level-set method for 3D process TCAD simulations on the feature scale. In the following, the two simulation tools are shortly introduced. Note that other tools exist as well, e.g., Sentaurus Process [134], but are not further discussed as they are out of scope.

Victory Process is a commercial process TCAD simulator developed by Silvaco [135]. The simulator allows to create a digital twin of electronic device manufacturing processes, such as etching, deposition, and oxidation. Victory Process uses a level-set engine based on hierarchical grids, allowing to represent different parts of material regions with varying spatial resolution. Details on the hierarchical grid data structure are presented in Chapter 2. Victory Process is written in C++ and pThreads and OpenMP are used for parallelization. The underlying simulation framework of Victory Process is the basis for the developed algorithms presented in Chapter 5, Chapter 6 and Chapter 7.

ViennaTS is a process TCAD simulator developed by the Institute for Microelectronics at the TU Wien, focusing on processing challenges for micro- and nanoelectronics [136]. It uses the level-set method on a hierachical run length encoding (HRLE) [99] data structure in combination with a sparse field approach. The open source topography simulator is written in C++ and uses OpenMP for parallelization.