(image) (image) [ Home ]

Emulation and Simulation of
Microelectronic Fabrication Processes

3.2 Chemical Surface Rate Modelling

The goal of this approach is to consider the physical processes which lead to the evolution of material interfaces. In order to obtain a surface velocity field \(\vec {V}(\vec {x})\), the chemical reactions on the surface have to be modelled in order to express their topographical effect on the material through etch or deposition rates \(R\). These reactions depend on the particles present at the considered surface site, expressed using surface coverages \(\Theta _x\) denoting the fraction of the surface site covered by a molecule of type \(x\). The distribution of coverages for each particle type is governed by the transport of this type in the gas phase through the reactor. Modelling this transport, the surface fluxes \(F_x\) for each particle type can be found in order to calculate the surface coverages \(\Theta _x\). Hence, as shown in Fig. 3.1, several different models have to be formed with the result of one model used as the input for the next. After all the physical models have been evaluated, the final result is the topographical change of the surface captured in the surface velocity field \(\vec {V}(\vec {x})\).

This section will discuss these models and their physical footing in detail. The models will be presented in the order in which they would be performed during the simulation of a semiconductor fabrication process.

(-tikz- diagram)

Figure 3.1: Flowchart of the chemical models used to generate the surface velocity field \(\vec {V}(\vec {x})\) and the information they provide for the next step (see arrows). First, the distribution of molecules impinging on the surface is modelled, which produces the incoming flux at every point on the surface. This flux is then used to model interactions with the surface in order to find how molecules are distributed leading to the surface coverage for each particle type. Using the number of molecules present at each surface site, their chemical reactions with the substrate are modelled, finally leading to the surface rates and thus the velocity field describing the topographical change. Applying these velocities to the initial surface using a single iterative advection step results in the final surface and the cycle starts anew.

3.2.1 Molecular Transport in Plasma Environments

The first step in the chemical description of a fabrication process, is calculating the rates at which different atoms, molecules, or ions impinge on the surface. In the context of process technology computer aided design (TCAD) modelling, these rates are referred to as particle flux and depend strongly on the different geometrical effects and transport phenomena inside the reactor [120, 121]. A thorough description of the thermodynamics of the fabrication process is necessary to capture all relevant phenomena. As already mentioned in Section 2.1.1, the reactor space is divided into a reactor scale and a feature scale region. Although reactor scale modelling is not the main scope of this work, the fundamental concepts are discussed in the following, focusing on the implications it has for feature scale modelling.

3.2.2 Reactor Scale Transport

As shown in Fig. 3.2, the simulation domain is split into reactor scale and feature scale regions, separated by the plane \(\mathcal {P}\), referred to as source plane. The reactor scale thus encompasses all space inside the reactor, whose dimensions are at the very least on the order of the wafer size, i.e. 300 mm for modern processes. A simple parameter used to determine the transport regime of particles is to consider their mean free path in an ideal gas:

\begin{equation} \bar {\lambda } = \frac {k_B T}{\sqrt {2} \pi d^2 P} \quad , \label {eq::meanFreePath} \end{equation}

where \(k_B\) is the Boltzmann constant, \(T\) is the temperature, \(P\) is the pressure inside the reactor and \(d\) is the collision diameter of the gas molecule [122], which for the most commonly encountered particles is around 0.4 nm [123].

(image)

Figure 3.2: Schematic representation of the traversal of neutral molecules and ions through the reactor and feature scale regions. While the motion through the reactor scale is dominated by random collisions with other molecules, the path through the feature scale region is dominated by ballistic transport. The directional distribution of neutral atoms and molecules \(\Gamma _{neutral}\) (blue) and ions \(\Gamma _{ion}\) (red) entering the feature scale region is shown as blue and red arrows, respectively. The molecular entities will then traverse the feature scale region in straight lines, only colliding with the surface.

In extreme cases, processes may operate at temperatures up to 1500 K and pressures down to 0.1 Pa [124], leading to a mean free path of \(\bar {\lambda } \approx \SI {290}{\milli \meter }\). For these types of processes, the exact reactor geometry and operation would have to be considered in order to generate accurate particle fluxes on the surface [125]. However, most microelectronic fabrication processes do not operate above 500 K or below 10 Pa, so the mean free path does not exceed \(\bar {\lambda } \approx \SI {10}{\milli \meter }\), meaning it is much smaller than the reactor region. Therefore, molecular transport is within the continuum regime and the directional distribution of particle motion is uniform, while the velocity follows the Maxwell-Boltzmann distribution [126]. Since particles are equally likely to move in any direction, the directional distribution of the particle flux crossing the source plane \(\mathcal {P}\) towards the wafer is described by a cosine [60]

\begin{equation} \Gamma _{neutral}(\theta ) = F_{neutral} \cos (\theta ) \quad , \label {eq::neutral_dist} \end{equation}

where \(\theta \) is the angle to the source plane normal and \(F_{neutral}\) is the neutral particle flux towards the wafer, usually found through experiments or reactor-scale simulations taking into account the properties of the specific reactor. In general, the flux \(\Gamma _x\) for any particle of type \(x\) may also include angular distributions for the energy of the particles or the location on the source plane. For neutral species, the kinetic energy is usually too low to significantly impact interactions with the surface. As this work focusses on structures well below the micrometre regime, it is also reasonable to assume that neutral flux distributions do not change significantly across the length scale of the simulation domain. Hence, the flux of neutral particles through the source plane can be described appropriately without considering the energy distribution of particles or the location on the source plane.

In numerous fabrication processes, biased plasmas [127] are used to accelerate ions towards the wafer and achieve focused particle trajectories. Essentially, the distribution of the ions becomes more directional as the plasma sheath potential increases and ion movement is not only governed by the random collisions with other particles, but also by the electromagnetic field in the reactor. Hence the particle flux for ions has a more focused directional distribution which can be approximated by considering the ion flux \(F_{ion}\) and a power cosine [128]:

\begin{equation} \Gamma _{ion}(\theta , \vec {x}_{\mathcal {P}}, E) = \frac {n + 1}{2 \pi } F_{ion}(\vec {x}_{\mathcal {P}}) \epsilon (E) \cos (\theta )^n \quad , \label {eq::ion_dist} \end{equation}

where \(n\) is the exponent of the cosine, \(\epsilon (E)\) is the ion energy distribution giving the fraction of all ions with energy \(E\) and \(\vec {x}_{\mathcal {P}}\) is the position on the source plane.

Depending on the type of plasma process to be simulated, different simplifying assumptions can be made. In a direct current (DC) plasma, if the thermal energy is small compared to the ion energy, all ions have the same energy:

\begin{equation} E = eV \quad , \end{equation}

where \(e\) is the elementary charge. Therefore, the ion energy distribution follows a Dirac delta function

\begin{equation} \epsilon (E) = \delta (E - eV) \end{equation}

and the exponent \(n\) can be approximated with [129]

\begin{equation} n = \frac {2E}{k_B T} \quad . \end{equation}

In focused ion beam (FIB) plasmas, it can be assumed that all ions have the same energy and direction and only depend on \(\vec {x}_{\mathcal {P}}\) as the ion flux is distributed according to a Gaussian [130] around the centre of the beam \(\vec {x}_0\)

\begin{equation} \Gamma _{ion}(\vec {x}_{\mathcal {P}}) \propto \exp \left ( {|| \vec {x}_{\mathcal {P}} - \vec {x}_0 ||_2^2} \right ) \delta (E - E_0) \quad , \end{equation}

where \(E_0\) is the ion energy of the beam.

The source distribution in molecular-beam epitaxy (MBE) reactors is more complex, as a ring of point sources is used to generate the source flux [131]. It is also possible that these point sources emit different materials, further complicating the arrangement [132]. It is important to consider the exact setup of the reactor in order to find the spatial and directional source flux for the feature scale which forms a small part of the wafer, since the directional distribution of the incoming flux changes drastically across the wafer [133].

Radio frequency (RF) plasmas have complex energy, spatial and directional dependencies and therefore plasma sheath simulations should be used to find an accurate expression for \(\Gamma _{ion}\) [134] in these types of plasmas.

3.2.3 Feature Scale Transport

By definition, the size of the feature scale is much smaller than the reactor scale region, namely on the length scale of the manufactured features, i.e. smaller than 10 µm. The mean free path of particles \(\bar {\lambda }\) defined in Eq. (3.1) is much larger than the size of the feature scale and thus the movement of particles is governed by ballistic transport, which can lead to shadowing due to the straight path of particles. Physically, each infinitesimally sized source plane element \(dA\) in Fig. 3.2 can be considered an individual particle source with properties governed by its source distributions \(\Gamma _x\). Due to the relatively slow movement of the wafer surface compared to the molecules in the gas phase, the surface can be treated as constant during the modelling of the particle transport. Practically, this means that it can be assumed that there is a constant flux from the source plane throughout the calculation of incoming particle fluxes on the wafer surface.

3.2.3.1 Boundary Conditions

Since it would be computationally unfeasible, and often impossible, to simulate the entire wafer, only a small portion of interest, such as a single feature, device, or circuit, is simulated. In order to still properly encompass the influence of neighbouring sections of the wafer, appropriate boundary conditions must be used for the simulation. Usually, structures would be repeated across the wafer, e.g. static random access memory (SRAM) cells in the memory region of a microprocessor. Therefore, when a particle leaves the simulation domain and hits a boundary, it must be treated accordingly. Using reflective boundary conditions, as shown in Fig. 3.3, the particle would simply be reflected off the domain boundary and therefore stay within the simulation domain. Periodic boundary conditions refer to the particle re-entering the domain through the other side, thus leading to the same description, as if the simulation domain was repeated periodically outside of the domain boundary.

3.2.3.2 Modelling Approaches

The two main approaches to modelling particle transport are numerical modelling using ray tracing methods and analytical calculations which require certain assumptions about the process to be made. In analytical models, certain properties of the initial geometry can be used to find a solution for the particle flux on the surface. One example for this approach is using Knudsen diffusion to model particle transport along straight trenches, which has been used successfully for specific processes [135, 136, 137]. However, deep knowledge of the modelled process is required and certain assumptions have to be made, so analytical models are not applicable to all processes and geometries. Furthermore, their general approaches may differ considerably, which is why they will not be discussed in detail here.

Ray tracing methods follow the ballistic trajectories of molecules in the feature scale region. However, due to the large number of molecules, not all can be simulated explicitly, but their trajectories are approximated using Monte Carlo (MC) methods. Each ray then represents potentially thousands of molecules, drastically reducing simulation runtime. Since the transport is simulated directly, the resulting distribution of molecules represents the expected physical result, given that the number of simulated particles is large enough. Since the transport simulation encompasses the trajectories of particles and the geometry directly, they are not limited to certain processes or initial geometries. Therefore, ray tracing is more appropriate for the modelling of general fabrication processes, where simplifying assumptions cannot be made easily.

The following sections will describe the numerical methods employed during ray tracing in detail.

3.2.3.3 Top-Down Flux Calculation

As suggested by its name, in the top down modelling approach, particles are represented using rays originating from each source plane element \(dA\). In the simplest case, all particles are adsorbed on the surface on impact and do not reflect or re-emit. Thus, the propagation of each ray ends somewhere on the surface. After all rays have been traced, the number of rays terminating at each surface element is counted, which is equivalent to the particle surface flux \(F(\vec {x}_s)\), as discussed in Fig. 3.1. Typical source plane fluxes for semiconductor manufacturing processes are on the order of 1016/(cm2 s) [110]. Therefore, in order to simulate a wafer area of 1 µm2 during one second of a fabrication process, roughly 108 rays must be traced. As most fabrication processes last substantially longer than one second, the large number of rays which would have to be traced makes it unfeasible to actually simulate each molecule individually. Hence, each ray usually represents thousands of molecules rather than just one. The particle surface flux \(F(\vec {x}_s)\) is then simply the product of the number of rays terminating at the surface element \(\vec {x}_s\) and the number of molecules per ray. Although representing multiple molecules with one ray leads to averaging of the surface fluxes, the result can be considered a good approximation as long as the total number of rays is much larger than the number of surface elements:

\begin{equation} N_{\text {rays}} \gg N_{\text {surface elements}} \quad . \end{equation}

To achieve a desired accuracy in the final result, the number of rays originating from the source plane should thus be related to the number of surface elements the rays may terminate on:

\begin{equation} N_{\text {rays}} \propto \frac {1}{\sigma _{\text {ray}}} N_{\text {surface elements}} \quad , \end{equation}

where \(\sigma _{\text {ray}}\) is the relative error compared to the analytical result for the surface flux, obtained by integrating over all source plane distributions \(\Gamma _x\). Intuitively, this can be verified by considering that, in order to achieve a numerical error of \(\sigma _{\text {ray}} \rightarrow 0\), an infinite amount of rays would have to be traced.

Given an expression for the source plane distribution \(\Gamma _x(\theta , \vec {x}_p, E)\), as in Eq. (3.2) and Eq. (3.3), MC methods are employed to generate the initial properties of each ray [138]. The starting position of each ray on the source plane can either be generated pseudo-randomly or chosen from a rectilinear grid of points on the source plane with a regular grid spacing. The latter is computationally more efficient and leads to comparable results, as long as the grid spacing on the source plane is comparable to the resolution of the surface the particles are incident on. Once all the necessary properties of the ray are known, its path to the surface can be traced using algorithms which have been extensively applied in computer graphics [139, 140, 141, 142]. The general idea of these ray tracing methods are shown in Fig. 3.3 for two particle source distributions, as they would be used in the simulation of a semiconductor fabrication process. The figure also highlights boundary conditions of the simulation domain, as well as possible reflections of molecules if they are not fully adsorbed by the surface.

(image)

Figure 3.3: Schematic depiction of rays being traced from the source plane to the surface using reflective boundary conditions. Each ray describes either neutral molecules (blue) or ions (red), governed by the source distributions \(\Gamma _{neutral}\) and \(\Gamma _{ion}\), respectively. Diffuse and specular reflections are shown for neutral species and ions, respectively.

The probability with which the particle will stick to the surface is referred to as the sticking probability \(S\) and may vary for each combination of molecule and material it is incident on. Since each ray represents numerous molecules, \(S\) can be used as a probabilistic factor rather than simply deciding whether the ray should be reflected or not. This means that the ray is always traced further as if the reflection occurred, with the flux carried by that particle reduced by a factor of \(1 - S\). When the flux reaches a pre-defined low threshold, the tracing of the particular ray is stopped. How particles are reflected or re-emitted from the surface depends on the specific chemistry and is governed by a reflection distribution \(\Gamma _{refl}\) describing the angle, flux and energy of the reflected ray, analogous to the source distribution \(\Gamma _x\).

Several algorithms have been proposed in order to perform ray tracing directly on implicit surfaces [143, 144, 145]. However, many of these algorithms stem from the field of computer graphics and are tailored towards visualisation applications. Hence, they are not necessarily applicable to the simulation of molecular transport. Furthermore, performing ray-surface intersection tests on implicit surfaces is computationally more expensive than on explicit surfaces [146], which is why molecular transport simulations are usually performed using explicit surface representations, such as triangulated surfaces [147]. However, the generation of a triangulated surface from a level set is expensive and thus undesirable. Hence, a different type of explicit surface should be used for the simulation of molecular transport using ray tracing.

In order for ray-tracing methods to work robustly, the explicit surface must be closed. There is no requirement for the mesh, however, to not contain self-intersections or other usually undesirable inconsistencies. Therefore, the explicit surface can be approximated using spheres of radius \(\Delta g\) centred at each grid point [148], or more accurately by discs forming tangential planes to the surface [60], as shown in Fig. 3.4. Spheres or discs can be extracted highly efficiently while still providing a robust, closed explicit interface for fast intersection tests of commonly used ray tracing algorithms. Since each sphere or disc is associated with a grid point, the results of the ray tracing step can easily be related back to the section of the surface described by the specific grid point.

(image)

(a)

(image)

(b)

Figure 3.4: Two ways to approximate an implicit surface efficiently by explicit shapes on active grid points (black), in order to simplify intersection tests for ray tracing. (a) Tangential line segments used to form an explicit approximation of the surface, as described in [60]; (b) surface approximated by explicit circles centred at active grid points, as described in [148]. The line segments and circles are replaced by discs and spheres in three dimensions.

In the top-down ray tracing approach, the physical behaviour of molecules moving through the feature scale is simulated and thus, all physical properties may be included in the simulation. Hence, it is also possible to include phenomena such as charged molecules in the feature scale being repelled by surface charges [149]. Furthermore, even complex reflective properties can be modelled using the top-down approach straight-forwardly. Specular reflections, for example, can be realised by considering the incoming angle of each ray with respect to the surface, which is found easily from the intersection test. Diffuse reflections employing any type of reflection distribution may also be used straight-forwardly, as discussed above. Crystal-direction dependent effects on the reflection or re-emission of molecules can also be included straight-forwardly by including the surface orientation at the intersection point [150]. Additionally, even particle-particle collisions can be modelled if the simulation domain is too large for the ballistic transport assumption to hold, such as high aspect ratio geometries like deep vias or trenches.

However, the inclusion of additional physical properties in the description of a process may result in a strong increase in the required computational effort. Although the top-down approach is often more computationally intensive than other forms of modelling molecular transport in the feature scale region, it is very closely linked to the physical behaviour inside the reaction chamber. Therefore, it is the most commonly encountered method in physical process simulations.

3.2.3.4 Bottom-Up Flux Calculation

The bottom-up approach to calculating the surface flux can be considered the reverse of top-down flux calculation, described in the previous section. Instead of tracing rays along the physical paths of the molecules, an element on the surface is considered and the paths of molecules are traced back to their origin. This is achieved by considering how much of the source plane is visible from a point \(\vec {x}_s\) [151, 152], as shown in Fig. 3.5. Numerically, the visible parts of the source plane \(\mathcal {P}\) are found by considering one surface element \(\vec {x}_s\) and iterating over all discretised points on the source plane. If the distribution is visible from \(\vec {x}_s\), its contribution towards the surface element is added to the total flux, according to each source flux distribution \(\Gamma _{src}\).

(image)

Figure 3.5: Schematic representation of the bottom-up flux calculation for modern transistor gate structures, using periodic boundary conditions, meaning the entire simulation domain is repeated at the boundaries. The black arrow indicates the direction used to find the direct flux incident on \(\vec {x}\) from a single particle source at \(\vec {x_\mathcal {P}}\) with a source distribution of \(\Gamma _{src}\). The flux of all visible source plane elements, indicated by the green arc, is summed to give the total direct flux on \(\vec {x}\).

The sum of all visible point fluxes results in the total molecular flux incident on the surface element at \(\vec {x}_s\)

\begin{equation} F_{\text {direct}}(\vec {x}_s) = \sum _{\vec {x}_\mathcal {P}} \Gamma _{src}(\vec {x}_s, \vec {x}_\mathcal {P}) Y(\vec {x}_s, \vec {x}_\mathcal {P}) \quad , \label {eq::direct_flux} \end{equation}

with the total flux being referred to as direct flux, since it does not include any reflections from other parts of the surface. The visibility function \(Y(\vec {x}_s, \vec {x}_\mathcal {P})\) is unity if the point \(\vec {x}_\mathcal {P}\) is visible from \(\vec {x}_s\) and zero otherwise. Since the exact flux emitted from \(\Gamma _{src}\) depends on the angle of emission, the distribution depends on the surface point and source plane point. This is indicated by the black arrow in Fig. 3.5. The angles at which source plane points are visible from \(\vec {x}_s\) are indicated by the green arc. All distributions on the source plane within the cone indicated by the dashed grey line thus contribute to the direct particle flux incident on the considered surface element.

In order to identify which points are visible and thus to find a description for \(Y\), adaptive visibility sampling can be employed. In this method, rays are launched towards the source plane in spatial directions spaced by pre-defined angles. If a ray hits the source plane, this part of \(\mathcal {P}\) is visible to the surface point. If the ray hits another part of the surface first, the source plane is not visible from \(\vec {x}_x\) in that specific direction.

The achieved resolution is given by the angle between two consecutive rays, since a transition from a visible to a shadowed section can only be resolved at this angle. The resolution can be increased adaptively, by launching additional rays in the direction of the observed transition [147]. Once a specific angular resolution has been achieved, all source distributions can be summed and the final flux can be calculated. If only a fixed resolution is required, numerous rays can be launched simultaneously, which can be performed efficiently on graphics processing units (GPUs) [153]. Therefore, sampling visibilities can be much more efficient than actually simulating the physical paths of molecules as in the top-down approach.

However, modelling reflections using the bottom-up method requires additional considerations, as they cannot be modelled directly from the source distributions. In order to find the reflected or re-emitted fluxes, an iterative method must be applied, which first solves for the direct flux \(F_{\text {direct}}\) which is then used to find the reflected and re-emitted fluxes. Each surface element is then considered as a particle source with distributions \(\Gamma _\text {refl}\) and \(\Gamma _\text {reem}\), as shown in Fig. 3.6. The source distributions include all the reflection properties based on the material and the type of the described molecule, such as energy, sticking probability, or angular dependence. Similarly to the direct flux defined in Eq. (3.10), the reflected and re-emitted fluxes incident on a surface element at \(\vec {x}\) are found by summing all contributions from other surface elements at \(\vec {x}’\).

(image)

Figure 3.6: Calculation of reflected or re-emitted fluxes using a bottom up technique with periodic simulation boundaries. The black arrow indicates the direction used to find the reflected and re-emitted flux incident on \(\vec {x}\) from a particle source at \(\vec {x’}\). The source distributions \(\Gamma _\text {refl}\) and \(\Gamma _\text {reem}\) define the flux emitted towards \(\vec {x}\). The total indirect flux at point \(\vec {x}\) is found by summing the flux from all visible surface points, highlighted by the blue arc.

For simplicity, the fluxes in the iterative method are named according to the current iteration, meaning the direct flux is labelled \(F_0\), the flux resulting from the first iteration of reflection and re-emission is labelled \(F_1\) and so forth. The i-th iteration of the secondary fluxes is given by

\begin{equation} \begin {split} F_{i+1}(\vec {x}) &= F_i(\vec {x}) - F_{\text {i,refl}}(\vec {x}) - F_{\text {i,reem}}(\vec {x}) \\ & + \sum _{\vec {x}’} Y(\vec {x}, \vec {x}’) \left [ \Gamma _\text {i,refl}(\vec {x}, \vec {x}’, F_i) + \Gamma _\text {i,reem}(\vec {x}, \vec {x}’, F_i) \right ] \quad , \end {split} \end{equation}

where \(F_{\text {i,refl}}(\vec {x})\) and \(F_{\text {i,reem}}(\vec {x})\) are the total outgoing fluxes from the surface element \(\vec {x}\). The arriving flux at the considered surface element \(\vec {x}\) is thus reduced by the total reflected or re-emitted flux and is increased by the arriving fluxes from other parts of the surface. This step is repeated until the desired accuracy has been reached or all particles have been absorbed, meaning that \(F_{\text {i,refl}}(\vec {x}), F_{\text {i,reem}}(\vec {x}) \rightarrow 0\) for all \(\vec {x}\). The number of iterations necessary to arrive at appropriate results strongly depends on \(\Gamma _\text {i,refl}\) and \(\Gamma _\text {i,reem}\). If they do result in a large portion of the arriving flux leaving the surface again, meaning that the sticking probability is low, more iterations will be necessary as more reflections or re-emission take place.

Since the flux at the surface element \(\vec {x}\) is preserved after each step of the iteration rather than individual particles and their properties, modelling specular reflections is not possible straight-forwardly. This type of reflections can only be approximated by assigning an average incoming and outgoing direction [154]. This leads to shortcomings when modelling ion-enhanced etching processes, which are heavily influenced by energetic ions and their reflections off the surface. Nevertheless, for processes with little influence from ions, the bottom-up approach can be used as a more efficient alternative to the top-down method.

3.2.4 Molecular Interactions with the Surface

The modelling of feature scale transport, discussed in the previous section, results in the incoming flux values on the surface for all molecular species involved in the fabrication process. The next step is to model how the particles which impinge on the surface interact with the substrate material to be adsorbed on the surface and thus influence the coverage of different molecular species on it [155]. If a material is covered by a species, it means that the specific molecules are adsorbed on the surface long enough to engage in chemical reactions leading to a topographical change on the surface. This is in contrast to reflected particles, which leave the surface immediately and are not able react with the substrate chemically.

The most general type of process, which encompasses physical and chemical surface interactions, is the ion-enhanced etching of a substrate. This type of process is used here as an example for the procedure of modelling interactions with the surface, as it may involve many different types of molecules and simultaneous deposition and etching. For simplicity, the particles which impact the surface may be grouped together according to the chemical effect they have on the substrate [156]. This results in four abstracted types of particles affecting the surface topography:

  • 1. Etchant: Leads to the removal of substrate material

  • 2. Passivation Species: Deposits on the surface to form a passivating layer

  • 3. Passivation Etchant: Removes passivating layer and possibly substrate

  • 4. Energetic Ions: Physically impacts the surface delivering kinetic energy

In order to appropriately model the process, the incoming flux of etchant \(F_e\), passivation \(F_p\), passivation etchant \(F_{pe}\) and ions \(F_i\) must be known to find the corresponding coverages \(\theta _e\), \(\theta _p\), \(\theta _{pe}\) [157].

The coverages \(\theta _e\) and \(\theta _p\) are normalised such that for each point on the surface, their sum cannot be greater than unity, while \(\theta _{pe}\) is defined as a fraction of \(\theta _p\). The coverages can therefore also be interpreted as the probability of the specific type of particle being present on the considered surface site. There is no coverage of ions, since ions do not deposit on the surface, but rather deliver kinetic energy for immediate physical reactions on the surface.

Since the transport of molecules is fast compared to the surface movement, it is assumed here that the coverages on the surface will reach a steady state immediately after the fabrication process is started [158, 159]. Based on this assumption, a linear equation can be set up for each particle type’s coverage [160], describing all relevant factors which contribute to the coverage.

Etchant Coverage

The etchant coverage is expressed as:

\begin{equation} \frac {d \theta _e}{dt} = F_e S_e(1- \theta _e - \theta _p) - k_{ie} F_i Y_{ie} \theta _e - k_{ev} F_{ev} \theta _e \approx 0 \quad , \label {eq::etch_coverage_rate} \end{equation}

where \(S_e\) is the sticking probability of etchant on the surface. \(k_{ie}\) and \(k_{ev}\) are the stoichiometric factors for ion-enhanced etching an chemical etching respectively, describing how much of one material, compared to its reactant, is needed in the chemical reaction. \(Y_{ie}\) is the ion-enhanced etching yield, which has a square root dependence on ion energy and incoming angle to the surface. It is usually modelled as [161]

\begin{equation} Y_{ie} = \left (\sqrt {E} - \sqrt {E_{th}} \right ) cos(\omega ) \quad , \end{equation}

where \(E_{th}\) is the threshold energy for ion-enhanced etching to occur and \(\omega \) is the angle between the incoming ion and the surface normal.

\(F_{ev}\) is the evaporation flux describing how much of the chemically etched material is transported from the surface through desorption and may be expressed as:

\begin{equation} F_{ev} = K_s F_e e^{- E_s / k_B T} \quad , \end{equation}

where \(K_s\) is a dimensionless reaction rate constant which defines how likely a reaction is to take place once the etchant comes into contact with the substrate and \(E_s\) is the binding energy between molecules of the substrate.

Passivation Coverage

The passivation coverage is given by:

\begin{equation} \frac {d\theta _p}{dt} = F_p S_p - F_i Y_p \theta _p \theta _{pe} - F_{depo} \approx 0 \quad , \label {eq::polymer_coverage_rate} \end{equation}

where \(S_p\) is the sticking probability of passivating molecules on the substrate and \(Y_p\) is the etching yield for ion-enhanced polymer etching. \(F_{depo}\) is the flux describing at which rate molecules are incorporated into the substrate during deposition in order to advance the surface.

Passivation Etchant Coverage

The coverage of passivation etchant is written as:

\begin{equation} \frac {d\theta _{pe}}{dt} = F_{pe} S_{pe} (1 - \theta _{pe}) - F_i Y_p \theta _{pe} \approx 0 \quad , \label {eq::polymeretch_coverage_rate} \end{equation}

where \(S_{pe}\) is the sticking probability of this type of particle on the passivation.

Physical considerations

The first terms in Eqs. (3.12), (3.15) and (3.16) describe the main mechanism with which the coverage is increased, namely through particle transport in the gas phase. Since only the molecules which actually stick to the surface without being reflected can engage in surface reactions, this term is also dependent on the respective sticking probabilities. The subsequent negative terms describe loss mechanisms which for etching occur only through the removal of the substrate through energetic ions or evaporation. In the case of deposition, the main loss mechanisms are ion-enhanced etching and the incorporation of molecules into the material through deposition.

Several additional phenomena which occur in specific processes could also be added. If surface diffusion of reactants play a role, a term describing this flux can be included. This flux is then simply obtained by considering the fluxes at close surface points and allowing for diffusive flow to the currently considered site.

In the simplest case, the particle fluxes do not depend on the coverages on the surface. In this case, ray tracing can be performed to find the fluxes which, in turn, are used to calculate the coverages for each particle type. However, in several processes, the reflective properties of particles may change drastically with changing surface coverages. This includes physical sputtering of the substrate through energetic ions, which is strongly dependent on etchant coverage and may be modelled as a type of re-emission. Sputtering, therefore, has to be considered during the modelling of molecular transport.

In order to resolve this cyclic dependency, an initialisation step is necessary for the coverages to settle into a steady state, as required for Eqs. (3.12), (3.15) and (3.16) to be valid. Once the coverages have converged to a value, the surface velocities can be calculated and the level set advanced. The coverages at the surface locations before the advection serve as a good approximation for the next particle transport modelling step [162], since the level set cannot be moved by large distances due to the CFL condition discussed in Section 2.4.2.2. Hence, this initialisation step is usually only necessary at the very beginnning of a simulation. However, it is crucial to perform this step for each different physical surface rate model applied.

3.2.5 Chemical Reactions of Molecules on the Surface

Once the fluxes and coverages on each surface point are known, they can be used to identify how the surface will evolve over time. In order to identify this evolution, the main processes which lead to a change in the surface must be identified. In the example of an ion-enhanced etching process, discussed in the previous section, there are several types of etching and deposition mechanisms which lead to topographical changes:

  • 1. Deposition of passivation

  • 2. Ion-enhanced etching of passivation

  • 3. Chemical etching of the substrate

  • 4. Ion-enhanced etching of the substrate

  • 5. Physical sputtering of the substrate

These reactions are depicted in Fig. 3.7, highlighting how the different fluxes and coverages discussed in the previous sections contribute to the topographical changes.

(image)

Figure 3.7: The five mechanisms leading to a topography change of the surface in the ion-enhanced etch process example. 1. Passivating molecules form polymers on sidewalls (green), which may be removed again through 2. ion-enhanced etching. 3. During chemical etching the substrate reacts with an etchant forming volatile etch products which thermally desorb from the surface through evaporation. 4. Ion-enhanced etching increases the chemical etch rate by breaking up bonds in the substrate, enhancing the formation of volatile etch products. 5. Energetic ions collide with the substrate and break the surface bonds to remove material physically without a chemical etchant through ion sputtering.

Passivation

The first two mechanisms only affect the passivation layer deposited in regions of low ion flux. This passivating layer therefore protects vertical features of the surface from being etched chemically. Eq. (3.15) can be rearranged to give the deposition flux used to advance the surface, which is simply the number of particles available for deposition on the substrate. The deposition rate of the passivation at a point \(\vec {x}_s\) on the surface is therefore given by

\begin{equation} v_{depo}(\vec {x}_s) = - \frac {F_{depo}}{\rho _p} = \frac {1}{\rho _p} \left ( F_p S_p - F_i Y_p \theta _{pe} \right ) \quad , \label {eq::surfaceModel_depoRate} \end{equation}

where \(\rho _p\) is the number density of the passivation material, meaning the number of particles per volume. The minus sign for the deposition rate stems from the convention of negative LS values denoting the inside of a material. The passivation coverage \(\theta _p\) is not included in the second term as deposition means that the surface is covered in passivating particles and therefore \(\theta _p = 1\). If there is a large flux of passivating particles at the surface, a passivation layer will be deposited, while a large ion flux results in the suppression of growth. If there is a previously deposited passivation layer beneath, the negative value means that this layer is etched. If the surface point \(\vec {x}_s\) is on the original substrate rather than on the passivation layer, then \(v_{depo}\) must not be negative, as passivation etchant cannot remove the substrate. In this case, \(v_{depo}\) has no physical meaning and can be ignored.

Substrate

The substrate can be modified through all of the above mechanisms with the exception of the ion-enhanced etching of the passivation layer. If deposition dominates and \(v_{depo} > 0\), no further considerations are necessary and \(v_{depo}\) can simply be applied at the specific point on the surface. If, however, \(v_{depo} < 0\) then etching of the substrate takes place. Each of the three etch mechanisms is expressed by one term in the final surface velocity:

\begin{equation} v(\vec {x}) = \frac {1}{\rho _{sub}} \left [ \underbrace {F_{ev} \theta _e}_\text {chemical etching} + \underbrace {F_i Y_{ie} \theta _e}_\text {ion-enhanced etching} + \underbrace {F_i Y_s (1 - \theta _e)}_\text {ion sputtering} \right ] \quad , \label {eq::substr_etch} \end{equation}

where \(\rho _{sub}\) is the number density of the substrate. Since all coverages affect \(\theta _e\), it must be found by solving Eqs. (3.12), (3.15) and (3.16). Chemical and ion-enhanced etching consume etchant and therefore result in a decrease of etchant coverage, while ion sputtering removes the substrate without affecting any of the coverages.

Encapsulating a large number of physical etch or deposition mechanisms, a single model could, in theory, describe a wide array of different fabrication processes. The choice of stoichiometric factors, threshold energies and other coefficients for each mechanism is crucial in order to properly capture the properties of a fabrication process. Since these physical constants entail all differentiating properties between two different processes, they must be considered with great rigour. Usually, these constants are found from experiments, reactor-scale or ab-initio simulations, or a combination of both. By adding more terms to the above model, any physical effect affecting surface movement could be included, although simple models, such as the one described here, are frequently enough to describe even complex processes [110]. It is therefore important to identify dominant etch or deposition mechanisms in order to ensure that the main properties of the fabrication process are modelled adequately.