Effect of Cu impurities on wet etching of Si(110): formation of trapezoidal hillocks

We simulate the formation of experimentally observed trapezoidal hillocks on etched Si(110) surfaces, describing their generic geometrical shape and analyzing the relative stability and/or reactivity of the key surface sites. In our model, the hillocks are stabilized by Cu impurities in the etchant adsorbing on the surface and acting as pinning agents. A model of random adsorptions will not result in hillock formation since a single impurity is easily removed from the surface. Instead a whole cluster of Cu atoms is needed as a mask to stabilize a hillock. Therefore we propose and analyze mechanisms that drive correlated adsorptions and lead to stable Cu clusters.


Introduction
Anisotropic wet chemical etching of silicon is a commonly used process for micromachinery fabrication. Using wet etching alone or combined with other methods one can make structures such as cantilevers [1], membranes [2], grooves [3], microfluidic channels [4], inertial sensors [5] and other micro-electromechanical systems (MEMS). As modern applications demand increasingly high precision at smaller scales, it is important to be able to produce defectless structures and high quality surfaces. In particular, the morphology and roughness of etched surfaces can affect device performance. As an example, the use of etching simulations in an injection mold application for the production of plastic razors using etched silicon as a molding surface has helped in understanding the step flow and step bunching origin of certain unwanted surface corrugations, thus allowing their elimination [6]. Similarly, the surface morphology developed for etched micro needle arrays can lead to larger friction and thus poorer detachment after drug delivery [7].
It is well known that in certain etching conditions hillocks may appear on etched surfaces [8]- [15]. Although the etching process should automatically remove any protrusions from the surface, a small mask can stabilize the needle-or pyramid-shaped structure of a hillock. It has been proposed that in order to have hillock formation during etching, four generic requirements should be fulfilled [16,17]: (i) the etch rate of the main orientation must be sufficiently high (the 'floor' must etch quickly); (ii) the etch rate of the facets of the hillock must be low (the hillock walls must be stable) [18]; (iii) the ridges of the hillock must etch slowly (the hillock ridges must be stable) [16]; and (iv) the apex of the hillock must be pinned (the apex must be stable) [19]. A masking agent is needed for the last condition to hold.
The spontaneous appearance of hillocks during etching means that unwanted, foreign agents are acting as masks. It is known that hydrogen released in the etching reaction can form bubbles on the surface and act as masks [12,20,21]. On the other hand, also impurities in the etchant may adsorb on the surface and act as pinning agents [13,22,23]. It has been 3 experimentally observed that Cu impurities in the etchant can lead to rough surfaces covered with hillocks [11,14,15] and first-principles calculations have verified that adsorption of Cu on a H-terminated Si surface is energetically favorable [24]. Other metal impurities do not have such a strong impact on the etched surface morphology [11,15]. As a result of the masking (or pinning) effect, pyramidal hillocks are observed on the silicon (100) surface [12,13,17,23] while on (110), zigzag profiles and trapezoidal hillocks may be seen [11,25,26]. The structure and formation of the zigzags and pyramidal hillocks has been analyzed and reproduced in simulations [13,23,27,28]. However, accurate analysis of the structure of the trapezoidal hillocks and the site-specific conditions for their stability, thus enabling their realistic simulation, has not been attempted before.
In this study, we examine the effect of Cu impurities on the morphology of Si(110) through simulations. We introduce a simulation method for silicon etching incorporating adsorbing impurities which may pin the Si atoms on the surface. We demonstrate, for the first time, how to produce trapezoidal hillocks in simulations and analyze the factors affecting their formation. The present study shows that there is a general model of hillock stability applicable to both the pyramidal and trapezoidal hillocks on (100) and (110), respectively, thus supporting the theory that the etching process and the surface morphologies of the two orientations are similar in spite of their different appearance, mainly resulting from differences in the underlying crystallographic structure, as described in [17].

First-principles calculations
We have calculated the adsorption energy of Cu on H-and OH-terminated Si surfaces at various surface sites using the atomic orbitals based SIESTA [29,30] implementation of density functional theory (DFT) [31,32]. The calculations were carried out using the PBE gradient corrected functional [33] and norm-conserving pseudopotentials. Basis set superposition errors were handled using counterpoise corrections [34]. The calculated systems were periodic twodimensional (2D) Si slabs with both sides terminated with H or OH to saturate dangling bonds. The Cu impurities are originally present as cations in solution. However, the copper ions are reduced near the surface and adsorb in the neutral state [11,24,35]. For a detailed description of the chosen calculation parameters and the used basis sets, consult [24]. We will also use some of the naming conventions from [24] in this study.

Kinetic Monte Carlo (KMC) simulations
We simulate the etching process using our own implementation of the KMC method [17,36]. The simulation of etching is based on an atomic description of the surface and the K-level search (KLS) [37] algorithm for choosing the Si atoms to be removed. The removal rates of individual Si atoms in the system are determined by their local neighborhood, i.e. the numbers of the first and second neighbor Si atoms in bulk and on surface [38]. This four-index classification makes it possible to distinguish and classify generally the different surface structures met in the simulations. Classification of the most important types is introduced in [39]. The removal rates used in this study follow closely the corresponding two-index values given in [17]. In addition to just simulating Si etching, we also want to model how Cu impurities in the etchant solution adsorb on the Si surface. The Cu simulation, described below, is joined to the Si etching simulation via a hybrid search tree as described in the appendix. The etchant is represented as a limited reservoir of Cu atoms which have certain probabilities to adsorb on the surface and then again desorb. In this study, we do not directly model the behavior of the Cu in the etchant nor do we account for surface diffusion. To incorporate impurities that may adsorb on the surface in the simulation, we must keep track of the constantly changing Si landscape and the potential adsorption sites on it. This is done by taking a 2D projection of the positions of H atoms and OH groups on the surface, thus generating a 2D point cloud, and constructing a Delaunay triangulation [40] for the points via the Bowyer-Watson method [41,42]. An example of such a triangulation is shown in figure 1. Unlike the example, however, the implemented triangulation has periodic boundaries. This triangulation is automatically updated during the simulation in order to monitor changes in the shape of the surface. Since a 2D projection is used, undercut geometries are handled incorrectly by the program. This is not a problem though, for we do not encounter such structures in this study.
The potential impurity adsorption sites in the simulation are defined to be the triangles in the surface triangulation. This serves two purposes: firstly, the surface triangulation automatically gives us a list of triangles which is equal to the list of adsorption sites. Secondly, the first-principles calculations have shown that the surface sites to which Cu can adsorb are usually located in between the H atoms on the surface [24] making the triangles a natural way to define the adsorption sites. The surface site types are distinguished according to the shapes of the corresponding triangles and each type is issued a rate for adsorption (and desorption) of Cu, where t i is the type of triangle i. It should be noted that some different sites, e.g. terrace sites C and D and kink/step sites F, F1, etc introduced in [24], are indistinguishable as triangles and have to be treated as a single type. Also, sites treated as quadrangles in the first-principles scheme are split into two triangles.
Once a Cu has adsorbed on the surface, it acts as a pinning agent. In the simulation this means that the removal rates of the Si which are connected to the Cu are multiplied by a pinning factor-10 −3 in our model. Here, 'being connected' means that a surface H atom bonded to the Si acts as a vertex in the triangle with the adsorbed Cu.

Nonuniform adsorption (NUA)
In our model, the role of Cu in hillock formation is to act as a pinning agent at the apex of the hillock. The pyramidal hillocks seen on (100) have very stable {111} facets and even a single impurity is enough to stabilize a hillock for a considerable time. The trapezoidal hillocks on (110) on the other hand have approximately {311} facets which etch away much faster than (111) even when they are stable compared to (110). This will lead to the masking agent being removed due to underetching much more quickly than in the case of the pyramidal hillocks on (100). Due to this mechanism, the stability of the facets also affects the stability of the apex. Furthermore, trapezoidal hillocks can remain stable for long periods of time only if the apex is masked by numerous Cu atoms or if there is a constant flow of new pinning agents to the apex. Due to simple diffusional principles, Cu clusters will not be stable if all surface sites have equal adsorption rates, as illustrated in figure 2(a): assume that a cluster of several Cu atoms is formed due to a random fluctuation and that it pins the apex of a hillock. The copper desorption and adsorption rates are proportional to the number of Cu atoms and the number of available surface sites, respectively. These are denoted by D and A, respectively, in figure 2, and the subscripts 'a' and 'f' refer to 'apex' and 'floor', respectively. Since the Cu density is at its highest at the cluster, we have D a /A a > D f /A f . Furthermore, if adsorption and desorption in the majority of the system (i.e. on the floor) are in balance, D f = A f and so D a > A a , i.e. there is a net flow of Cu atoms away from the Cu cluster and the mask dissolves. Therefore, in order to have stable masks of numerous atoms, the adsorption cannot be uniform. Instead, there must be a higher adsorption rate at the vicinity of the Cu cluster than elsewhere in the system, as shown in figure 2(b). If the adsorption rate is enhanced enough, A a D a and the net flow of Cu to the cluster will be positive (or zero at the saturation size of the Cu cluster) and a stable mask is formed.
One could suspect that the nonuniformity in adsorption rates is simply due to the high number of possible adsorption sites [24]. However, the surface sites found at the apices of hillocks are abundant also on the rough (110) surface, i.e. on the floor. Enhancing the adsorption to some particular surface site type will not only increase the adsorption rate at the apex (A a ) but on the floor as well (A f ). As a result, the etch rate of the floor is reduced and the hillock will not grow. Thus, an additional mechanism must make the Cu atoms prefer adsorbing near each other and we refer to these as NUA mechanisms in this paper. Below, we describe three schemes for realizing such behavior in simulations and comment the physical relevance of the models.

Interaction enhanced adsorption (IEA).
An attractive interaction between Cu impurities may result in NUA. If it is energetically favorable for two Cu atoms to sit next to each other on the Si surface, the Cu may cluster and thus form a stable mask which stabilizes the apex of a growing hillock.
To test the feasibility of this idea, we use DFT to calculate the energy of a (221) oriented supercell with two Cu atoms on distant (111) terrace sites (C and D in [24], see figure 3(a)) and compare it with one where the Cu atoms are on adjacent but otherwise similar sites ( Figure 3(b)). The energy of the paired configuration is found to be 0.9 eV below the one with two distant Cu atoms, showing a strong attraction. Also, the calculated charge densities of the two configurations are clearly different as the two Cu atoms bond when brought close to each other. Thus, clustering due to Cu-Cu interactions is possible.
We call our implementation of this mechanism in KMC the IEA. In this scheme, we simply monitor all surface triangles and keep a record of how many of their neighbors hold an adsorbed Cu atom. Here, we define neighboring triangles to be such that they share a common vertex. Based on the number of neighboring Cu, n i , we multiply the adsorption rates in (1) by an enhancement factor f IEA For simplicity, we mostly use a constant factor. The IEA is a completely local mechanism for NUA as it only acts according to the immediate neighborhood of each site.

Height enhanced adsorption (HEA).
While IEA makes the impurities find each other, NUA may also be achieved by making the impurities prefer adsorption at the apices of hillocks. In practice, this would mean monitoring fluctuations in the surface morphology and amplifying them by correlating the impurity adsorptions to peaks in the Si landscape. There is more than one feature of the surface that may be used for identifying the apices, but the simplest one to use is probably the height of the surface. The scheme using this approach is named HEA. Similarly to IEA and (2), we define an enhancement factor f HEA . However, now the factor is a function of surface height, h i , measured from the lowest point of the surface: In order to ignore very small features of the surface and to also prevent the enhancement from becoming too large, we use a truncated linear function for f HEA , Unlike IEA, HEA is a completely global scheme since it depends on the absolute height. We implement HEA out of theoretical interest, as the basic concept of the method is different to IEA. HEA effectively enhances the adsorption of copper at the highest locations of the surface. Initially, a fluctuation in the height profile of the surface leads to a first maximum and adsorption is enhanced on it, effectively pinning the surface at that stochastically chosen location for the rest of the simulation. Physically, HEA may be viewed as a way of mimicking the effect of a boundary layer in the etchant. However, such layers may be very thick in a real etchant and thus the associated effect should be weak. This renders the physical relevance of the effect questionable. Therefore, we are mostly interested in qualitatively describing the hillock formation resulting from this type of mechanism.

Activity enhanced adsorption (AEA).
In order to introduce a method which is not completely local, like IEA, nor global, like HEA, we implement the AEA scheme. The motivation for using such a method comes from [28], where activity monitoring was successfully applied in the simulation of diffusion effects in the study of step bunching during etching. AEA is a phenomenological method in which we keep track of the regions where Cu has adsorbed and adjust the adsorption rates to make the adsorptions correlated. By increasing the adsorption probability at the neighboring regions of the sites where Cu adsorption has already occurred, Cu clustering is made possible stochastically without actually introducing an attractive force between the coppers. The method mimics the action of an attractive force whose magnitude would grow with the cluster size without actually simulating the force directly. The principle is similar to that of HEA with the exception that AEA monitors and amplifies fluctuations in Cu density while HEA is coupled to fluctuations in surface height. The technical implementation is as follows. We define an (abstract) activity field, A i , over the surface. Every time a Cu adsorbs, the activity at the neighbors of the adsorption site is increased. On the other hand, there is a rate at which the activity is being constantly lowered. The adsorption rates are determined by the activity in a similar way as in IEA and HEA, Here, we use a linear factor The activity is limited both below (must be non-negative) and above (can be A max at maximum). The purpose of A max is to limit the cluster sizes and prevent the attraction between Cu atoms and Cu clusters from strengthening perpetually over time. Note that if the activity was decreased at each desorption instead of gradually, the method would be equivalent with IEA. Figure 4(a) is an SEM image of a trapezoidal hillock seen in experiments as presented in [11]. The characteristic shape of these hillocks is clearly seen in the image: the quadrangular base is approximately diamond-shaped with the diagonals parallel to the axes of symmetry of the Si(110) surface. The hypothesized structure of the trapezoidal hillock is shown in figure 4(b). Thinking in terms of {111} layers, the observed shape consists of stacks of such layers shaped as triangles. The facets are approximately {311} oriented [11,17] and consist ideally of dihydride   [11] in (a) should ideally consist of stacks of triangular (111) layers [17], as in (b). This hypothetical structure is realized in the simulations, as in (c) (using IEA). Trapezoidal hillocks have not been previously realized in simulations due to the complexity of the specific details required for the fulfillment of the four general conditions for stable hillock growth (section 1). First of all, the hillocks can form only if the etch rate of the (110) floor is considerably higher than that of the {311} facets ((1) and (2) in figure 4(d)). In KOH etching this means that a high KOH concentration is required. On the atomistic level relevant in a simulation, the relative etch rate of monohydride steps found on the floor must be higher than that of dihydride steps found on the facets. More difficult to overcome in simulations is the fact that the edges and the apex of the hillock must also be stable for the hillock to survive ((3) and (4) in figure 4(d)). Achieving stability of the apex requires masking by, e.g. impurities. How this may be realized is discussed in section 2.3. In principle, impurity particles can also accumulate on the facet edges stabilizing them [13]. However, the simulation models applied in this study do not strongly drive such behavior and therefore we fit the Si etch rates to make the edges always relatively stable. (This can be thought of as manually pinning specific sites.) It turns out that only the removal rates of two special Si types appearing at the edges need to be specifically defined. These types are the previously unclassified (2, 1, 0, 5) and (2,1,3,3) in the four-index scheme from [39]. We lower their removal rates from the corresponding two-index values given in [17] by a factor of 10 −3 .

Structure of trapezoidal hillocks
Once the four aspects are taken into account, the trapezoidal hillocks are finally seen in simulations, as shown in figure 4(c). The {111} layers are clearly visible in the figure and the simulations confirm that the assumed structure in figure 4(b) is indeed correct. The precise shape, size and density of hillocks depend on the etching conditions used in the simulation and the applied scheme of NUA, as we discuss in the following sections.

Cu adsorption rates
3.2.1. Ideal rates. First-principles calculations have shown that it is energetically favorable for Cu to adsorb on the Si surface with the adsorption energy being approximately proportional to the number of bonds the Cu can form [24]. Due to the energy hierarchy, there are a few sites where adsorption is especially favorable. Using the naming convention of [24], the most favorable sites are the F, H and I as well as A if OH groups are present on the surface. These special sites are shown in figure 5. The F and H sites appear mostly on horizontal dihydride steps and also at some kinks, A appears on monohydride steps and I appears at kinks.
To study how pinning at these special sites affects the morphology, we run a series of calculations using ideal adsorption rates. That is, we systematically go through all subsets of {A, F, H, I} and assign high adsorption rates for the chosen sites while all other sites have a very low adsorption rate. This way practically all available Cu in the system will be adsorbed at sites of chosen types. The desorption rates are all set low and desorption occurs through the underetching route (the appendix), i.e. when the neighboring silicons are removed. The etching conditions are adjusted so that hillocks form when adsorption is allowed on all the sites A, F, H and I. (In this case it means a reservoir of about 200 Cu atoms for IEA and AEA and 1000 atoms for HEA while the 75 × 75 nm 2 surfaces hold 150-200 thousand sites.) We were unable to produce hillocks for any combination of adsorption rates (not limited to the sites discussed here) when no NUA scheme was applied. Table 1 shows which combinations of adsorption sites result in stable hillocks. We see that hillocks are easiest to create within AEA, as only the combinations where both A and H sites are denied adsorption end up without hillocks. For HEA, also the case where adsorption on only A sites is allowed produces a flat surface. IEA is the most demanding scheme. Using it, hillock growth is achieved only if adsorption on site A as well as F or H is allowed. This delicateness in IEA is due to locality: the interactions only play a role if there are two surface sites with a reasonable adsorption probability next to each other. If we prevent adsorptions at key sites, this seldom happens.
We conclude that pinning of certain key sites is sufficient for the trapezoidal hillocks to grow. More specifically, all of the NUA methods result in stable hillocks if adsorption is allowed on all of the sites A, F, H and I. The AEA scheme produces hillocks even if adsorption is denied from all other sites except A or H. This suggests that hillocks induced by Cu impurities are due to pinning of Si atoms at these key sites. There is a wide range of adsorption rates for which hillock growth is observed, and therefore knowledge of the precise adsorption rates at the various surface sites is not required to qualitatively simulate the formation of the hillocks.

Rates estimated from first-principles results.
In order to link our KMC simulations with the DFT data, we want to estimate Cu adsorption rates from first-principles results. Thorough calculation of adsorption rates for all possible sites is too expensive, though, and we must resort to simplifications. This suffices, since in the previous section we showed that the exact values of the impurity adsorption rates are not needed in order to qualitatively simulate the process of hillock formation.  [24] with H-termination (squares) or OH-groups present (circles). X marks uncategorized sites. The dashed line shows a rough split to potential adsorption sites and relatively much more inert sites.
In order to determine adsorption rates for KMC, we need activation energies for adsorption. The activation energies, E, allow us to approximate the adsorption and desorption rates, k, within transition state theory as where k B is the Boltzmann constant and T is the temperature. The prefactor ν is left as a parameter in order to synchronize the timescales of Si etching and Cu adsorption simulations, since our Si removal rates are only relative, not absolute. We assume one ν ads to hold for all adsorption rates and one ν des for all desorption rates. As examples, we calculate the adsorption energy barriers for sites D and F, located on a (111) terrace and a dihydride step, respectively, as shown in figure 5. This is done with DFT using the drag method. The transition state energies found for these sites are 1.6 and 1.5 eV above the energies of the adsorbed states, respectively. These values are directly the activation energies for desorption. The corresponding barriers for adsorption, given by subtracting the adsorption energies from the aforementioned numbers, are 0.4 and 0.1 eV, respectively. In order to estimate the adsorption rates of other surface sites, we make the crude assumption that the transition state energies for all sites are 1.5 eV above the adsorbed state. In other words, we assume that adsorption occurs most rapidly on the sites that have the highest adsorption energies.
The resulting activation energy hierarchy is shown in figure 6. We distinguish between purely H-terminated surfaces (squares) and sites where at least one OH-group is present as a vertex of a surface triangle (circles). For cases with OH where DFT results are unavailable (open circles), we estimate the energy to equal that of an E or F site, whichever is closer to the energy of the same site with H-termination. This is done in agreement with the observation that the presence of OH-groups makes the adsorption energy hierarchy more uniform [24]. Since the adsorption rates depend on the activation energy exponentially, the estimated hierarchy essentially allows adsorption at sites F, H and I as well as A if there is OH-coverage and denies it elsewhere. (This is the same group of sites examined using ideal rates in the previous section.) We use this set of rates in the following simulations.

Influence of the NUA mode on hillock formation
Having shown that trapezoidal hillocks can form in simulations, we proceed to analyze the formation process and resulting morphologies. Based on the simulations with ideal adsorption rates, we know that the implemented NUA schemes lead to qualitatively different results. Furthermore, the methods include parameters which affect the morphologies. Therefore, we analyze the different cases separately.

IEA.
In IEA, as far as hillock formation is concerned, the decisive parameter is the interaction strength. Here, we use the enhancing factor f IEA = exp(E/k B T ), where the parameter E is effectively the decrease in activation energy due to the interaction. (However, here E is not limited by the actual activation energies.) In figure 7, we show surfaces simulated under IEA with two different values for E (0.5 and 0.8 eV). The figure also shows the adsorbed Cu atoms as white circles. We see that for a 0.5 eV interaction, the hillocks are unstable even though there are masking Cu clusters present. This is because the masks themselves are unstable. When one Cu atom adsorbs, the interaction is strong enough to make most of the other Cu impurities cluster around it. Some Cu atoms are bound to adsorb elsewhere though, and clusters start to form also around them. As the number of Cu atoms is limited, there is a competition between the clusters. The new clusters are growing on the (110) floor where there are plenty of favorable adsorption sites available while hillocks have already started growing beneath the older clusters reducing the number of good adsorption sites. Especially, there are fewer pairs of adjacent, favorable adsorption sites on the apices of  (3)) is varied as shown.
hillocks. Because of this, it is most likely that the Cu atoms adsorb to the newest clusters and the older ones eventually dissolve as the Cu atoms desorb.
When we strengthen the interaction to 0.8 eV, a trapezoidal hillock forms. Now the interaction is strong enough to pull all available Cu in one cluster and it will take longer until any competing clusters appear. The interaction is also strong enough to make adsorptions happen at other sites besides just the most favored ones. Thus, the mask and the hillock are stable. If we increase the amount of Cu, in the case of a weak interaction the masks remain unstable. For strong interactions, the resulting Cu clusters are larger, but we still only see one. This is an effect of the applied model: all the Cu atoms in the etchant can 'see' the entire simulated surface at all times. Therefore a strong interaction attracts all Cu atoms in the system to the same cluster. In reality, it will take some time before dissolved coppers find a cluster on the surface and during this time they may adsorb somewhere else and act as seeds for new clusters.

HEA.
In the height-dependent scheme, adsorption is not emphasized only at the very edge of a Cu mask on top of a hillock, but all over the hillock instead. This means that the impurities will be scattered on the hillock with highest density on the apex. In order to get sufficient masking at the top, much more impurities are needed than in IEA.
The efficiency of the HEA scheme depends mostly on the coefficient a in (3), which determines the strength of the nonuniformity in adsorption. For a very low a, no hillocks can be seen. When a is somewhat increased (to around 5 nm −1 ), hillocks appear but only in simulations with a very high Cu concentration since only a small fraction of the impurities find the hillock apices. The resulting hillocks are also quite round as there is Cu scattered all over the surface. When a is further increased, the Cu atoms seek the tops more efficiently and lower Cu concentrations are sufficient for producing the hillocks. The results of simulations with both a high and a low a at a medium Cu concentration are shown in figure 8. We also see in the figure that impurities seek not only the top of the hillock, but also the edges to some extent, forming a cross-shaped pattern. This is due to the edges having more favorable adsorption sites than the facets.  (5)) are varied as shown.
Since the height-dependent method uses a global height scale, the highest hillock, which was the first one to appear on the surface, is the most efficient in attracting Cu atoms and so it is also the most stable. This way, although several hillocks can form, one or two dominate the others and grow to be the largest hillocks. The simulation with a high a in figure 8 is an example of this. When the maximum height h max in (4) is reached, the whole part of the hillock above this level will have the same enhancement factor making the adsorption essentially uniform for that region. So, the peak of the hillock above h max becomes unstable and a plateau forms. Thus the hillocks cannot grow higher than h max .

AEA.
The AEA method is a compromise between the local interaction driven model and the global height-dependent scheme. Although each adsorption amplifies the activity field in an area which is equal to the interaction range in IEA, the activity can remain after the impurity has desorbed. Thus there is no requirement to have the most favorable adsorption sites next to each other like in IEA. Instead, it is enough that the favorable sites appear frequently in the same region. This makes Cu masks more stable. Importantly, several Cu clusters can coexist so we see many hillocks in the simulations. On the other hand, the method is not absolute like HEA, so none of the hillocks dominates the system. Figure 9 shows the simulated morphologies for two enhancement factors b in (5) and two frequencies of activity truncation. For the weaker enhancement factor (top row in figure 9) and the short truncation interval (left column) the result is quite dispersed adsorption on the surface. If truncating is done less frequently, Cu forms clusters more effectively and some hillocks are seen. This happens because when truncation is performed frequently, the activity around new clusters stays relatively low for a long time and the nonuniformity in the adsorptions is weak. If truncating is less frequent, only few adsorptions are needed to reach a high activity and cluster formation proceeds. Adjusting the enhancement factor higher has a similar effect. It makes Cu clusters grow faster due to faster rise in adsorption rates and so allows the hillocks to start growing earlier. This is seen in the lower row (high b) of figure 9 having larger hillocks after similar etch time than the top row (low b). The number of hillocks increases as either parameter is increased due to the number of non-clustered Cu atoms decreasing. Indeed, AEA appears to be the most flexible of the NUA methods since not only does it produce hillocks with the widest range of adsorption rates it also allows the hillock density to be adjusted by changing the simulation conditions.

Discussion
Of the proposed NUA methods, the IEA mimics the attractive interaction seen in first-principles calculations. Simulations with this scheme demonstrate that stable Cu clusters are a necessity in order to produce trapezoidal hillocks. On the other hand, due to the particular characteristics of the model, IEA simulations never show multiple hillocks at the same time. They typically produce a single cluster, which sometimes becomes unstable against the formation of new clusters at the flatter regions on the floor surface. These features have motivated us to introduce the height and activity-dependent methods, HEA and AEA, respectively. All three methods can be used for simulating cluster growth. HEA typically requires many impurities since it is not very efficient for this purpose. In comparison, AEA provides a balance between the number of impurities and the typical size of the clusters (which actually remain stable) distributing them homogeneously and stochastically over the surface. This activity-based scheme has proven especially useful as a general implementation of fluctuation amplifying feedback. Besides mimicking attraction, the method is also known to be efficient in the simulation of diffusion effects [28]. By comparing the three methods, we conclude that the most essential requirement for the formation of trapezoidal hillocks is the onset of copper clustering.
Experimentally, it is known that the trapezoidal hillocks appear once the Cu concentration exceeds a certain threshold level [11]. This feature is seen also in the simulations although the threshold depends on the applied method. HEA requires a higher Cu concentration before hillocks appear compared to the other methods. IEA works even for very low Cu amounts. In fact, at high concentrations, IEA generates large irregularly shaped Cu clusters and the resulting hillocks are also quite irregular in shape. (This is due to the limited size of the simulated systems. Were the hillocks as large as those seen in experiments, the irregularities resulting from the shapes of the masks would be insignificant.) AEA works for a wide range of Cu concentrations as the amount of Cu mainly affects the number of hillocks in the system.
The facets of the experimentally seen trapezoidal hillocks have been reported to be {311} planes [11]. However, this is an approximation. The exact orientation is determined by the relative stability of the surface Si atoms at various parts of the hillock. It is especially important that since the symmetry of the hillock is two-fold instead of four-fold like in the pyramidal hillocks seen on (100), there are two different types of facet edges (see figures 4(b) and (d)). The relative stability of these edges affects the shape of the trapezoidal base of the hillocks. This is in agreement with the experiments where the trapezoidal base and the facets of the hillocks can be stretched and/or rounded [35]. The stability analysis of the hillocks also verifies that the conditions of hillock formation proposed in [17] and discussed in section 1 as well as figure 4(d) hold for trapezoidal hillocks: the facets and facet edges must be stable with respect to the (110) floor, since trapezoidal hillocks form only for combinations of Si removal rates for which this is the case. (In simulated etching of Si(110) where the edges of a potential trapezoidal hillock, shown in figure 4(d) (3), are unstable, it is still possible to realize rectangular hillocks with {111} and {100} facets as described in figure 30(a) of [17]. This further demonstrates that the geometry of the trapezoidal hillocks is determined by the kinetics of the etching process.) Furthermore, the need for Cu clusters shows that the apices of the hillocks must be stabilized.
In KOH etching, the requirement that the facets of a trapezoidal hillock should be stable compared to the (110) floor is fulfilled only at a high KOH concentration, in agreement with the experiments where trapezoidal hillocks are observed only at these high concentrations (section 3.1) [17]. Additionally, this study has shown that adsorption of Cu on the surface site A, appearing predominantly on monohydride steps, plays an important role in stabilizing the apices of trapezoidal hillocks on (110) (section 3.2.1). According to the DFT calculations, adsorption of copper on site A only occurs at partially OH-terminated sites, as demonstrated by figure 6. Therefore, Cu atoms will adsorb on type A sites only at high KOH concentrations. Thus, both the stability of the facets as well as the pinning of the apex require a high KOH concentration. Note that the pyramidal hillocks on (100) only appear at low KOH concentrations since the (100) floor etches rapidly only in such conditions [17]. However, since coppers may adsorb on the sites F, H, and I even at a low KOH concentration, they can mask the apices of the pyramidal hillocks by pinning these sites.
Our model assumes that the masking effect of Cu impurities is due to the impurity atoms blocking the Si at the adsorption sites they occupy. However, it is possible that the effect is due to a more complex process. It is known that hydrogen bubbles on the Si surface can also act as micromasks. If adsorbed Cu atoms affect the wetting properties of the surface, the lifetime and size of the bubbles is also affected [21]. This way Cu could induce pinning indirectly. If this is the case, experiments could reveal a change in the density and size of H bubbles as a function of Cu concentration. It has also been proposed that hillocks form during etching because of pinning by semipermeable particles [13,23]. The hypothesis was that there are particles which adsorb, reduce etching rates and stay on the surface even as etching proceeds. The simulated Cu clusters act in this manner. The semipermeable nature of the cluster stems from the fact that individual Cu atoms can be removed, but the other remaining impurities attract new Cu atoms and so the 'particle' stays on the surface despite underetching.

Conclusions
We have simulated etching of Si(110) using a KMC scheme including etchant impurities that may adsorb on the surface and act as pinning agents. In this study, we focus on the role of Cu as the pinning impurity. For the first time, the analysis of the structure of the trapezoidal hillocks and the site-specific conditions for their stability has led to their simulation. The abstracted geometrical structure of the simulated hillocks agrees completely with the theoretical analysis and closely resembles the experimentally observed morphologies. We observe that clustering of Cu is a necessary requirement for the formation of the hillocks. If clustering occurs, the sitespecific adsorption of Cu at the most favorable sites resulting from DFT calculations is sufficient to drive the mechanism. First-principles calculations suggest that there is an attraction between the Cu atoms. Monte Carlo simulations incorporating such an attraction (IEA) show that it can indeed lead to the desired Cu clusters. Besides an interaction model, we have also implemented  Flowchart presentation of the used hybrid search tree. The choices between Si and Cu events are made stochastically according to the relative probabilities of the events. If we remove a Si atom to which a Cu is connected, the Cu desorbs automatically due to underetching. and analyzed two other computational schemes which drive Cu clustering and result in hillock growth. The present study shows that there is a general model of hillock stability applicable to both the pyramidal and trapezoidal hillocks on (100) and (110), respectively. adsorption. The total rate for Si removal is the sum of removal rates of all surface Si. Here, the total rate for Cu adsorption is the average of the adsorption rates of all unoccupied surface sites multiplied by the amount of Cu in the etchant-reservoir. The total desorption rate is the sum of the desorption rates of all adsorbed Cu atoms on the surface. Individual KLS trees are kept for surface Si and adsorption sites, i.e. unoccupied surface triangles. Both trees are updated after each physical event. In addition to thermal desorption, adsorbed coppers are also removed if a silicon to which the Cu atom is connected is removed.