Surface morphology during anisotropic wet chemical etching of crystalline silicon

The rich variety of micron-scale features observed in the orientation-dependent surface morphology of crystalline silicon during anisotropic wet chemical etching is shown to have its origin at the atomistic scale. Realistic Monte Carlo simulations show that the pyramidal hillocks on Si(100) are the result of local stabilization of distributed apex atoms by (metal) impurities from solution. In the absence of this stabilization, shallow round pits are formed on Si(100) as a result of the anisotropy between (one layer deep) pit nucleation and (isotropic) step propagation. It is also concluded that nosed zigzag structures at vicinal (110) are the combined result of misaligment and the etching anisotropy, showing that the nucleating mechanisms of morphologically related structures such as pyramidal hillocks and nosed zigzags are not necessarily the same. The simulations confirm that the formation of (one layer deep) triangular and hexagonal pits on exact Si(111) and of polygonal (saw-shaped) and straight terraces in vicinal Si(111) depends on the relative rate of [12̄1] and [1̄21̄] step propagation and on the misorientation of the surface with respect to Si(111).


Introduction
Anisotropic wet chemical etching remains the most widely used processing technique in silicon technology. Used in combination with a multitude of other processes, it has a wide range of applications in microelectromechanical systems (MEMS), including pressure [1,2], acceleration [3], angular rate [2] and gas-flow sensors [2,4], actuators [5], nanoprobes [6,7], nanowires [7], micromirrors [8], laser cavities [9], optical switches [10,11], aligment grooves [12] and microvalves [13], to mention only a few. Its wide presence is not only due to its ease of use and low cost, but also due to the fact that it provides fairly smooth surfaces with no physical damage to the bulk. However, in the past few years, the roughness of the etched surface has begun to play an increasingly important role as the size of the devices is continuously reduced. The performance of a growing number of micromechanical devices requires very smooth surfaces and the precise production conditions need to be found. It has become apparent that there is a need to improve the understanding of the process in general and, in particular, of the mechanisms that lead to the characteristic morphological features of the etched surfaces.
During recent years, much effort has been dedicated to the characterization and understanding of the surface morphology during wet chemical etching, both experimentally [12], [14]- [27] and theoretically [17], [28]- [34]. Specific studies of the most frequent surface inhomogeneities, such as pyramidal hillocks [12,17], [20]- [23] (figure 1(a)) and shallow round pits [24,25] (figure 1(b)) on Si(100), nosed zigzag structures [28] (figure 1(c)) on (vicinal) Si(110) and triangular pits [19,26,27], [30]- [34] (figure 1(d)) on Si(111), have provided important insights into the problem. As an example, there exists extended consensus that the atomistic nucleation of pits on Si(111) [31,34] is enhanced by the existence of linear dislocations [26] accumulated in the surface region during thermal oxidation [27] and that the development of a triangular or hexagonal pit geometry depends on the relative values of the microscopic reaction rates of the atoms sitting at the [121] and [121] step families [27,31]. Recently, a mechanism explaining the dependence of these shapes on concentration/coverage has been suggested [30]. Similarly, there is also agreement that the pyramidal hillocks on Si(100) are necessarily nucleated by a stabilizing mechanism, although the nature of the mechanism itself is still in debate.
Systematic studies of the dependence of the surface morphology on orientation [14]- [16] have provided a more general perspective on the problem, establishing morphological relations between some of the surface features such as the pyramidal hillocks on Si(100) and the nosed zigzags on Si(110) or the formation of polygonal-stepped (figure 1(e)) and straight terraces (figure 1(f)) in vicinal Si(111). Also, from a more general perspective, the specific features of the surface morphology in anisotropic wet chemical etching of silicon are only an example of a wider range of morphological features observed in other etched or grown materials. An improved understanding of the mechanisms generating these features in wet chemical etching may help in other disciplines.
The purpose of the present work is to stress the fact that most of the morphologic features in wet chemical etching are directly the result of atomistic processes. This is an important observation that has been traditionally overlooked, especially in the case of the formation of pyramidal hillocks, whose origin has been related to larger objects such as re-deposited islands, silica particles or hydrogen bubbles, as discussed in section 3.5. We will show that a realistic model based on simple atomistic considerations can generate simultaneously (with the same parameters) pyramidal hillocks on Si(100), nosed zigzags on (vicinal) Si(110), polygonal 100. 4 and flat terraces on vicinal Si(111) and triangular pits on exact (111). In doing this, we will learn about the specific amplification mechanisms that provide the link between the atomistic processes and the macroscopic morphologic features. In the particular case of the pyramidal hillocks and round pits on Si(100), the temperature and concentration/coverage dependence will be examined. Based on the fact that pyramidal hillocks and shallow round pits on Si(100) are observed under different experimental conditions, we will use two sets of parameters. Interestingly, it will be shown that the general features of the surface morphology at the other orientations do not depend critically on the choice between these two.
As the paper deals with all the relevant morphologic features that are observed in wet chemical etching, we have considered it most suitable to postpone the presentation of the available experimental/theoretical information about each specific surface feature to the beginning of the corresponding section where that particular feature is considered. The paper is organized as follows: in section 2 we describe the atomistic model that is used to perform the simulations. In section 3 we consider the morphology of Si(100). We first review previous experimental studies on pyramidal hillocks in section 3.1, present our results in sections 3.2-3.4 and discuss them in relation to other proposed mechanisms in section 3.5. We similarly consider the formation of shallow round pits in sections 3.6 and 3.7. The morphology of the other orientations is considered in section 4, including exact and vicinal (110) (section 4.1) and (111) (section 4.2). We present our conclusions in section 5.

Atomistic model
Anisotropic wet chemical etching is a non-equilibrium process in which both the microscopic roughness and morphology and the macroscopic orientation-dependent etch rate are determined by the relative values of the microscopic (atomistic) reaction rates. Gosálvez et al [35] have shown that the origin of the (large) differences in site-specific rates is found in two microscopic mechanisms: the weakening of backbonds following OH termination of surface atoms and the existence of significant interaction between the terminating species (H/OH). The weakening of the backbonds depends only on the total number of hydroxyls attached to the two atoms sharing the bond and is independent of the particular distribution of the OH groups between the two atoms [35]. The energy of a bond between an atom terminated by i OH groups and an atom terminated by j groups (i, j = 0, 1, 2, 3) can be written as where 0 ≈ 2.7 eV is the bond energy between two bulk atoms and ≈ 0.4 eV is the energy reduction for every OH group that is attached to either atom. Correspondingly, the total bonding energy for a surface atom with n first neighbours is simply the sum of the energies of the n bonds: where we have considered the most general case, in which the target atom is terminated by m OH groups (m ≤ 4 − n) and the j th first neighbour ( j = 1, 2, . . . , n), having itself n j first neighbours, is terminated by m j OH groups (m j ≤ 4 − n j ). The other microscopic mechanism of major importance in wet chemical etching, namely, the interaction between the surface terminating groups (H/OH), occurs only in the presence 100.5 of indirect second neighbours [30,36]. Due to these interactions, hydroxyl termination of the target atom (and its first neighbours) involves additional energy terms, not taken into account in equation (2). As a result, the total (local) energy of a surface atom can be expressed as the sum of three terms [36]: where E bonds is the energy of equation (2) and (e TA OH/H + e TA OH/OH ) ( (e FN OH/H + e FN OH/OH )) symbolically denotes the total energy from the interactions between the OH groups terminating the target atom TA (the first neighbours FN) and H and/or OH terminating the indirect second neighbours of the target atom TA (first neighbours FN). The geometrical restrictions to hydroxyl termination in the presence of indirect second neighbours is a manifestation of the important role of steric hindrance in anisotropic wet chemical etching. In the present model, the source of steric hindrance is identified as the (H/OH-terminated) indirect second neighbours.
Note that, although the parameters and 0 used for describing the bonding energy are fixed by the first-principles ab initio study [35], the interaction energies e TA,FN OH/OH and e TA,FN OH/H can be used as tunable parameters in order to describe different etchants. Once an etchant is chosen, its concentration is described in the model by the amount of surface coverage by OH groups, denoted θ . Following [36], we choose e TA OH/H = 0.2 eV, e FN OH/H = 0.0 eV, e TA OH/OH = 0.9 eV and e FN OH/OH = 0.7 eV in the present study. Similarly, hydroxyl termination of the target atom and its first neighbours (i.e. OH clustering) is assumed for all surface coverage values [36]. For the second neighbours of the target atom, one type of termination or the other will randomly occur depending on the OH coverage value θ .
The dynamics of the surface consists in random removals of surface atoms with probabilities where the activation energy E α is defined as Here, p 0α and E α are parameters describing the different surface atom types (α = 1, 2A, 2B, 2C, 3A, 3B). Note that the local energy E is calculated using the same expression (equation (3)) for all site types independently of the value of α. The function max(0, E − E α ) is used to conform with the Metropolis algorithm [37]. Following the discussion in [36] and the notation used in surface studies of Si(111) [34], we consider the following surface site types.
• Type 1, singly bonded atoms: trihydrides (TRI); also referred to as kinks.  Note that the atoms of type 0 are included for completeness since they can occasionally appear in connection to the formation of overhangs. This is, however, a rare event in the simulations and has no measurable effect on the evolution of the surface. These atoms are removed (with unit probability) as soon as they are encountered. Thus, in this model the surface contains six types of atom. Note also that due to the different possible combinations of the terminating species H and OH around a surface site, the activation energies E α = E − E α will take different values even for the case of atoms of the same type α.
The six pairs of parameters ( p 0α , E α ) for types 1, 2A, . . . , 3B can be determined from comparison to experiment. The idea is to choose the parameters so that the relative values of the etch rates of a number of surface orientations (six, in principle) agree with those from an experiment. By adjusting the parameters p 0α , the simulated etch rates will shift up/down in an Arrhenius plot. Similarly, the slopes of the etch rates can be controlled by tuning the parameters E α . Alternatively, it is also possible to choose the parameters ( p 0α , E α ) based on comparison of the simulated surface morphology with that from experiments. An example of this approach is provided by the present study. A similar approach is used in [31]. Note that since different sets of parameters ( p 0α , E α ) lead to different Arrhenius dependences for a set of surface orientations, these parameters can be used to represent different etchants.
Based on the fact that the pyramidal hillocks and the shallow round pits on Si(100) are observed under different experimental conditions, we will consider two sets of parameters, as indicated in table 1. Case A leads to the formation of hillocks and case B to the formation of round pits. As we will see, both cases provide the formation of nosed zigzag structures on Si(110) and triangular pits on Si(111).

Morphology of (100)
The morphology of Si(100) after wet chemical etching is characterized by the appearance of (i) near-pyramidal hillocks [12,17], [20]- [23] and (ii) shallow round pits [24,25]. Typically, hillocks and pits are reported to appear under different experimental conditions, not together. In both cases the density of surface defects (either hillocks or pits) increases as the concentration of the etchant is reduced [12,38] and the density of hillocks increases with temperature [12]. As we shall see, pyramidal hillocks are the macroscopic result of an atomistic process in which surface atoms are stabilized by (metal) impurities from the solution against removal by the etchant. The hillock nucleation process is independent of the way in which the surface as a whole is etched away (e.g. pit nucleation followed by step propagation, atomistically rough etching etc). Shallow round pits appear when the stabilization mechanism is not present and as a result of the anisotropy of the etching process between downward propagation of the surface (pit nucleation) and isotropic propagation of steps parallel to the surface.

Pyramidal hillocks
The formation of hillocks has been reported for a number of etchants such as hydrazine [39], ethylenediamine pyrocatechol (EDP) [40,41], ammonium hydroxide (NH 4 OH) [42], tetramethyl ammonium hydroxide (TMAH) [22] and, especially, potassium hydroxide (KOH) [17,18,20,21,28], [43]- [53]. At low magnification, the hillocks appear as regular pentahedra composed of four lateral {111} crystallographic planes lying on the (100) base plane. Under closer examination, the hillocks reveal a more complicated morphology, departing from the basic pyramidal shape in a way which depends on the particular experimental conditions. For instance, near-pyramidal hillocks after KOH etching were described by Tan et al [23] as being bounded by eight {567}-family planes, a few degrees off from the four {111} planes. The near-pyramidal morphology can be considered as the result of {111}-terrace bunching due to the relatively fast removal of material from the top and the four near-110 edges as compared to the slow removal of material from the {111} planes themselves (section 3.2).
At any given time, the hillock sizes are inhomogeneous, indicating that these surface protrusions are created continuously through the etching process [12,54]. The times required to grow hillocks to a size of several micrometres are rather long, in the range of 1-10 h, depending on etching conditions [17,20,47]. For conditions producing hillocks, if the exposure time is long enough, the surface is eventually completely covered with pyramids (complete surface texturization) [12]. At low concentrations, complete texturization of the surface can be achieved in relatively short times (tens of minutes) [54]. Incomplete texturization in the steady state is also possible for experimental conditions that set a limit to the maximum possible pyramidal size [18].
The surface density of pyramids decreases with increasing etchant concentration [12,18,40,43,47] and increases with increasing temperature according to an Arrhenius behaviour [12]. To our knowledge, the present work offers the first consistent explanation of this behaviour with concentration and temperature.
The formation of pyramidal hillocks cannot be correlated to the use of masks since hillocks are reported both on blank wafers and inside the open regions of masked wafers [18,22]. They cannot be eliminated by careful surface preparation (including cleaning of surface contaminants) prior to etching or by vigorous stirring during the etching process [23]. For conditions producing hillocks, stirring reduces the hillock density [23] and the hillock size [22]. Stirring also affects the uniform distribution of hillocks over the surface so that less repeatable patterns are obtained than without stirring [22]. The use of magnetic stirring also has a beneficial effect [48]. Similarly, Baum et al [20] report that the application of ultrasound reduces the average size of the pyramids from 2-3 to about 0.1 µm in 2 M KOH at 60 • C.
According to Yang etch experiments [17], the pyramidal hillocks are neither caused nor stabilized by linear or planar defects in the bulk of silicon. This is compatible with the fact that dislocations are known to accelerate the local etch rate in Si(111) rather than to decelerate it [26]. Pyramid formation has been correlated to both water quality and etchant purity [48]. In cases where extremely poor quality water is used, pyramid-covered surfaces are obtained and, in some 100.8 instances, minimal etching of the surface is observed. Small concentrations of metal impurities in solution (ppm or even ppb levels) are known to induce pyramid formation [45,55,56]. Similarly, the concentration of silicon in the solution is also reported to induce hillocks [46].
By re-etching in the same or different experimental conditions, the size of the near-pyramidal hillocks can be reduced and, in some instances, the hillocks can be completely removed [18,23]. During the re-etch process the pyramidal/near-pyramidal hillock morphology can change in some cases from having a near-square bottom to an octagonal base [51], the protruding region being enclosed by eight {331} facets [23]. The transformation of the near-{111} facets into {331} planes (as well as other non-octogonal transformations [52]) occurs by lateral attack of the near-110 edges.
For any etchant, high concentration may be used to avoid hillock formation. However, this will typically reduce the anisotropy between {100}, {110} and {111}, and high concentrations are typically only used to remove already existing hillocks [12,18,40]. Bressers et al [50] and Xia and Kelly [51] report that pyramids can be suppressed by the addition of oxidizing agents (such as ferricyanide, chromate, permanganate or oxygen) or by application of a potential positive with respect to the open circuit potential (OCP). Similarly, Campbell et al [48] report the absence of hillocks on the addition of isopropanol (IPA) to the etchant. The main effect of the oxidizing agents and the potential is to replace the Si-H site expected normally in chemical etching by an Si-OH site through chemical or electrochemical oxidation of a surface intermediate of the etching process. According to Haiss et al [53], in addition to Si-H, the predominant surface intermediate formed during etching contains a deprotonated hydroxyl group-O − (as in Si-O − ). Electrochemical indications of the important role of the OH − ions for the removal of hillocks were reported earlier on by Tan et al [12].

Simulation results for pyramidal hillocks
In this section we show the results of simulations using the parameters of case A (table 1) at different temperatures and coverages. Figure 2 shows the typical time evolution of the (100)surface morphology. Initially, the morphology is characterized by nucleation and growth of pyramidal hillocks. At any time, pyramids of different sizes are observed, reflecting the fact that hillocks are nucleated continuously during the etching process [52]. For long etching times, complete surface texturization by the hillocks is obtained, in agreement with experiment [12,54]. Note that, as in the experiments, the simulated hillocks are not geometrical objects with flat facets joining at a single apex. Instead, they contain rather irregular edges along the 110 -directions (figure 2(g)) and closer inspection shows that several points simultaneously share the function of the apex (figure 2(h)). Also, the flat (111) facets have a complex terraced structure, as shown in figure 2(h). This feature corresponds well with the fact that, after pit nucleation, the {111} surfaces are etched by a step propagation mechanism resulting in the formation of terraces, as will be shown in section 4.2.
A series of snapshots is used in figure 3 to show the typical nucleation mechanism for the formation of a pyramidal hillock. The hillocks are nucleated as a result of the stabilization of the initial apex atom due to the attachment of two metal impurity atom/ions from the solution. Note that, although not shown, the attachment of the two impurities is sequential (i.e. it does not usually occur simultaneously). In the simulations, we assume that the impurities interact with the singly coordinated silicon atoms, replacing them instantaneously. Effectively, the metal impurities are represented by the singly bonded silicons. The number of metal impurities on the surface is controlled by adjusting the parameters E 1 and p 1 entering the removal probability of singly coordinated silicons. This approach for the incorporation of metal impurities on the surface will prove satisfactory for the description of the temperature and coverage dependence of the surface density of hillocks (section 3.4). Note that the microscopic process effectively controlling the nucleation of a pyramid is the disclosure of singly coordinated atoms. In principle, these disclosures can occur in various multiple ways through the removal of two or three-bonded atoms and, effectively, they cannot be associated with a single atomistic process.
The removal of the stabilizing impurities will result in the destabilization of the apex atom and, thus, of the hillock. However, the hillock may be restabilized again by the new attachment of impurities to one or more surface atoms at the top region. This explains the appearance of distributed apex structures, in agreement with experiments [21]. In this way, if there are enough impurities and their removal is not very frequent, the pyramids can grow from their bottom down as the wafer surface evolves downwards at a faster rate than the pyramid tops. Note that, in addition to stabilized tops and fast removal of bottom-surface atoms, the formation of a pyramid requires two extra conditions: the stability of {111}-facet atoms (TM or type 3A) and the stability of ridge atoms (SM or type 3B). This is automatically provided in our model by the backbond weakening effect of the terminating OH groups. The atoms of type 3A are backbonded to three bulk atoms through strong bonds since, at most, the bonds are weakened by one OH group. In turn, the atoms of type 3B are bonded to one bulk atom and two surface atoms so that, typically, two of the bonds are weakened by two OH groups each and the third bond is weakened by only one group. Thus, under all possible conditions, this makes the removal of ridge atoms more probable than that of facet atoms. This explains the appearance of terraced structures at the interfacet region of the pyramids and the formation of macrosteps (in the sense of more than one layer thick) as observed in figure 2(h). This also explains why, in the absence of stabilizing impurities, the pyramidal hillocks are removed by direct attack to the ridge regions, as figure 4 shows. In this figure, a pyramidal hillock grown in similar conditions as those for figure 2 is re-etched at higher coverage and lower (non-zero) impurity concentration, keeping the temperature constant.

Macroscopic etch rate
The fact that the surface is eventually texturized (i.e. covered by hillocks) has a dramatic effect on the value of the macroscopic etch rate, as figure 5 shows. The etch rate experiences a large change (of about two orders of magnitude) as a function of time. This has major implications for the experimental measurement of the etch rate, since the usual determination as the ratio of the total etched depth to the total time becomes an ill defined quantity. This ratio will be physically meaningful only for two (extreme) cases: (i) when the surface texturization is still far from being complete, in which case it will provide a good estimate of the initial larger slope, and (ii) when the total time considered is so long that, in relative terms, complete texturization was reached soon after the initiation of the etching process. In this case, the ratio will provide an approximation to the smaller slope.
For intermediate etching times, the estimation of the smaller slope by means of a straight line joining the origin and the end point of the curve (i.e. the previous ratio) will provide a poor approximation. Note that, since the complete texturization of the surface may take hours in real experiments [17,20,47], it will be difficult to know whether the previous ratio can be used as an estimate of the etch rate unless the surface is monitored to detect the onset of complete texturization. In the worst case, the use of such etch rates in Arrhenius diagrams may lead to the erroneous determination of the activation energy. The existence of a gradual change from short-to long-term roughness and etch rate is reported experimentally [47].

Dependence on temperature and coverage
The density of hillocks after etching during a fixed period of time increases strongly with increasing temperature (figures 6(a)-(d)) and decreases with increasing coverage (figures 6(e)-(h)), in agreement with experiments [12,18,40,43,47]. In order to obtain a quantitative measure of the temperature dependence of the number of hillocks independently of the etch duration, we consider the rate of hillock formation, defined as the slope of the linear increase in the number of pyramids as a function of time before complete surface texturization ( figure 7(a)). Note that the density of hillocks should increase linearly with time, reflecting the fact that, for a fixed temperature and concentration, hillock nucleation is a completely random process and, as such, the longer one waits the larger the number of observed processes. As shown in figure 7(b), the rate of hillock formation follows an Arrhenius dependence on temperature, in agreement with experiment [12]. Note, however, that the hillock densities produced in the simulations are large (by orders of magnitude) compared to the experimental densities reported. This is due to the fact that, from a computational point of view, the realization of a statistically accurate determination of the dependence of the hillock density on temperature and coverage becomes too expensive for very large systems containing smaller densities of larger pyramids. At first glance, the fact that the density of hillocks increases with temperature seems to contradict the nucleation mechanism (figure 3). Since the removal probability of the stabilizing impurities increases as the temperature is raised, one would expect the number of stabilizing events to drop and, accordingly, a reduction in the number of hillocks. However, raising the temperature will also cause an increase in the removal probabilities of the other species as well. Thus, taking into account that it is possible to conceive of a situation in which the number of disclosures is increased with temperature more strongly than the number of impurity removals, while not seriously compromising the stability of the hillocks. This precisely describes the present case. Note that, according to this process, the concentration of impurities in solution should decrease with increasing temperature as a result of their consumption in the termination of the increasing number of apex-becoming atoms. Note also that the density of pyramids will decrease with increasing temperature (as reported in [18]) if the number of disclosures increases with temperature more slowly than the number of impurity removals. According to the previous discussion, the value of the activation energy for the rate of hillock formation ( figure 7(b)) should be a complicated average of both the microscopic processes that lead to singly coordinated atom disclosures and the removal of impurities. In principle, the relative importance of the different microscopic processes can be evaluated following a similar method to that used in the identification of the contributions to the activation energy of the etch rate [57]. The fact that the activation energy of the rate of hillock formation is similar for the two coverage values considered ( figure 7(b)) is explained in terms of the reduction experienced by the removal probabilities of all particle types with increasing coverage. An example of this effect is shown in figure 7(c), where the etch rate of the surface reduces with increasing coverage without significant change in activation energy.

Comparison to other nucleation mechanisms
A number of mechanisms for the formation of pyramidal hillocks have been suggested in the literature. We briefly compare the most significant cases with the present approach. One common feature of all proposed mechanisms (including the present one) is the recognition that pyramid formation can only be due to the stabilization of part of the surface by the absorption or deposition of material which can behave as a micromask. However, the nature of the masking agent and the way how it gets to the surface constitutes the difference between the different mechanisms and is a matter of debate.
(a) Hydrogen bubbles. As the etching process is accompanied by vigorous formation of hydrogen bubbles on the evolving surface [20,21], [48]- [51], the growing bubbles have been suggested to act as micromasks before detachment, inducing the formation of the pyramidal hillocks [20,21,48,49]. In principle, this idea is supported by the fact that pyramid-free surfaces are obtained in oxygen-saturated etchants [48], a phenomenon attributed to the rapid reaction of oxygen with hydrogen, resulting in reduced hydrogen bubble formation. The opposite occurs in solutions saturated with nitrogen, where pyramid formation is increased. In addition, the formation of circular arrangements of hillocks on vertically etched wafers [18] serves as a record of bubble formation, growth and movement on the surface, strongly supporting bubble involvement in hillock formation. In our opinion, the hydrogen bubbles may only occasionally nucleate the hillocks. The reduction of the hillock density in oxygen-saturated etchants [48] is not due to a reduction in the number or size of the hydrogen bubbles, but to an increase in the amount of surface coverage by OH groups, as concluded in [51] for other oxidizing agents and strongly suggested by the formation of a deprotonated hydroxyl reaction intermediate Si-O − in [53]. The increse in OH coverage results in fewer pyramidal hillocks, as shown in section 3.4. This is the same effect as obtained with an increase in the concentration of the etchant or with the use of oxidizing agents and biasing potentials [50,51]. We recognize that the H 2 bubbles have an important role in the amplification of the hillock size and in the generation of circular arrangements of hillocks on the surface, as strongly evidenced in [18]. Since the bubbles are mobile, they can act in coalition with metal impurities as momentary stabilizing agents, amplifying the effects of the impurities. The fact that the circular arrangements of hillocks in [18] are accompanied by the formation of random hillock distributions on the regions not affected by the bubbles supports our argument that the bubbles can modify and amplify the effects of more microscopic nucleation mechanisms. (b) Polymerized residues/reaction pr oducts/Si O 2 precipitates. Seidel et al [58] reported the observation of polymerized residues on the surface of silicon and suggested their relation to the formation of hillocks. However, Tan et al [12] find no evidence of residue localization at the top of hillocks. These residues are found on areas with no hillocks and even under etch conditions where no hillocks form at all. Similarly, no systematic relation between bulk SiO 2 precipitates and hillocks is found by Schröder et al [18], although occasional stabilization by SiO 2 micromasks cannot be completely ruled out. (c) Re-growth mechanism. Tan et al [12] suggest a re-deposition (or re-growth) mechanism based on under/super-saturation of the reaction products locally at the surface. Supersaturation leads to re-deposition at a local site if the transfer of the reaction products away from the solid-liquid interface is slow compared with its production. Occasionally, redeposition of Si clusters creates islands which then nucleate the hillocks. Further growth of the hillocks is not caused by the growth reaction (i.e. the reverse to the etching reaction), but by etch anisotropy. This mechanism is supported by the fact that increased concentration of silicon in the solution is reported to induce hillocks [46]. In our opinion, since the etching process is far from equilibrium, the reverse reaction is very unlike and can only explain a small fraction of the formed hillocks. In addition, the assumption of intrinsic re-growth due to super-saturation cannot be reconciled with the fact that external factors, such as metal impurities, are known to have an effect on the morphology, the roughness and the etch rate (see (e)). Should the impurities act as catalysts of re-growth, then this process should rather be referred to as assisted re-growth. (d) Semipermeable silicate particles. Nijdam et al [17] and van Veenendaal et al [28] elaborate further the nature of the stabilizing agent and its interaction with the surface. They suggest the hypothesis of semipermeable silicate particles that are attracted to the pyramid ridges. The most intriguing feature of these particles is their permeability: they adhere to the surface 100.16 but simultaneously allow the atoms directly underneath to be etched away, although at a lower rate than the atoms not covered by these particles. The hypothesis that the protruding ridges act as sinks of tiny particles would explain the experimental result that the mesoscopic etch rate of a pyramid edge is lower than the experimental etch rate of a Si{110} surface. In our opinion, the gathering of silicate particles at pyramid ridges might work in a two-dimensional world, where the pyramid ridges become the pyramid 'faces'. However, in the three-dimensional case, additional motion of the silicate particles parallel to the surface would be required. In addition, the explanation of how the particle provides etching of material underneath and simultaneously adheres to the surface (to ensure perpendicular propagation) in terms of the semipermeability property, would imply mesoscopic porosity of the silicate, in order to allow for diffusion of hydroxyl groups and/or water molecules through the particle (to produce attack) and diffusion of the products back to the solution. Thus, the particles should be rather large and, therefore, easily observable in experiments, which is not the case. (e) Impurity micromasks. Landsberger et al [22] and Campbell et al [48] have pointed out that precursors of micromasks are probably present as impurities in the etching bath. In the experiments of Campbell et al, pyramid formation is correlated to both water quality and etchant purity. However, they do not elaborate further on this issue, focusing instead on the hydrogen bubble mechanism (see (a)). More recently, Tanaka et al [55,56] have shown that very small concentrations of metal impurities such as Cu and Pb have clear effects on the etch rate and the morphology of (100) and (110). In addition, the simulations of van Veenendal et al [28] show that pyramidal hillocks are formed when the apex atom is protected by a two-atom micromask. By combining these experimental and computational results with our simulations, we conclude that the atomistic metal impurities are the principal cause for the formation of pyramids. Note that oxide micromasks and/or re-growth of silicon can occasionally provide the necessary nucleation mechanism and that the hydrogen bubbles have an important role, providing an effective amplifying mechanism. The technical question of whether the impurity atoms replace the singly coordinated silicons, simply terminate them or interact in some more complicated manner needs still to be addressed. Also, the possibility that the impurities may catalyse the re-deposition of silicon should be considered.

Shallow round pits
In the absence of pyramidal hillocks, the morphology of the (100) surface is typically built up from shallow round pits [14,25]. Sometimes the surface has a very smooth appearance to the naked eye, but under closer inspection, the surface reveals a very rough nature at the micron scale [38]. The formation of shallow round pits takes place during a timescale of hours [25]. It is well known, although not understood, that the round pits can be avoided by increasing the concentration of the etchant [38], in which case macroscopically smooth surfaces are observed. By producing surface indentations with different indent shapes and indent orientations prior to etching, Shikida et al [24] have shown that the round shape of the pits is inherent to the etching process. The indents always evolve towards a round shape independently of the initial indent shape and orientation. No systematic experimental studies on the dependence of the density and/or size of the round pits on temperature have been found. The formation of round pits is sometimes thought to be associated with hydrogen bubbles, especially due to the geometrical match of the two objects. However, hydrogen bubbles can only 100.17 retard etching locally, not accelerate it (section 3.5(a)). The round pits have also been related to bulk microdefects [25], although no conclusive evidence is provided and the origin is still debated. We will present here computational evidence that no external agents are required for the formation of shallow round pits. These surface inhomogeneities are the result of the anisotropy of the etching process between pit nucleation and step propagation, for the particular case in which step propagation is isotropic. The formation of round pits or macroscopically smooth surfaces depends on the amount of anisotropy existing between pit nucleation and step propagation.

Simulation results for shallow round pits
In this section we show the results of simulations using the parameters of case B (table 1) at different temperatures and coverages. Figure 8 shows three snapshots of the typical time evolution of the surface morphology. At the initial stage, the time evolution is characterized by nucleation and growth of shallow pits (figure 8 (top)). As the pits grow larger, eventually interacting with each other (figure 8 (centre)), new pits are initiated inside them, providing a mechanism to recycle the process (figure 8 (bottom)). As the average size that the pits may reach before interacting depends on the pit density, we consider larger densities of smaller pits in our simulations as compared to experiments [25,14] as similarly done in the case of the pyramidal hillocks (section 3.2). As the surface profiles show, the depth of the simulated pits is about 1% of the diameter. In comparison to the experimental etch pits, whose depth is 5-10% of the diameter, the simulated pits are too shallow. This is a result of the choice made to obtain large pit densities in small systems.
The mechanism of pit formation is a combination of (one layer deep) pit nucleation and isotropic step propagation, as shown in figure 9. Initially, one layer deep pits are nucleated randomly over the surface due to the removal of type 2A atoms (TD). The planar size of the nucleated pits grows due to step propagation, a process with a lower activation energy than the nucleation of the pit. Note that at any time two or more pits may meet, leading to the formation of a non-stable larger pit which becomes rounded soon after. From time to time, one (or more) pits may be nucleated inside an already existing pit, close to the expanding front i.e. away from the pit centre. As a result, the pit becomes deeper and locally more robust against interaction with neighbouring pits, therefore preserving its shape around that region in the collisions. As more nucleation events occur close to the perimeter, the pit becomes more robust against collision with other pits, becoming larger and deeper, but always shallow. Also, as the pit grows in diameter, the probability that new smaller pits are initiated inside it is increased, this leading to a selffeeding process in which the larger the pit becomes, the larger it will still grow. When two larger pits with similar size and depth meet, they combine, producing non-stable non-round pits which eventually become round. Note that the proposed mechanism explains the experimental fact that larger pits contain smaller countless copies inside them [28]. Figure 10 shows that the density of pits increases with temperature, as expected. According to the mechanism of formation (section 3.7), the number of nucleated pits is proportional to the removal of type 2A atoms (TD). Accordingly, the density of round pits should follow an Arrhenius behaviour with an apparent activation energy related to the activation energy of these atoms. We have not pursued this issue any further in this study.     Figure 11 presents the results of etching at higher coverage the same initial systems as those of figure 10 for the same temperatures and times. The fact that the pits can be removed by increasing the concentration of the etchant is in agreement with experiments. The reason for this behaviour is that isotropic step propagation at low coverage becomes anisotropic at high coverage, as shown by the formation of one layer deep elongated pits in figure 10(a). This occurs due to the fact that, as the coverage is increased, the removal of type 2C sites (HSD) becomes more disfavoured than the removal of type 2B atoms (VSD), introducing an anisotropy that is not present at low coverage. This anisotropy in step propagation at high coverage is the result of one extra indirect second neighbour for HSD, as was discussed in [36]. Note that the removal of type 2A sites (TD), which nucleate the pits, occurs far less frequently than the removal of VSD.

Temperature and coverage dependence
For smaller increases in coverage ( figure 11(b)), the surface becomes atomistically rough, showing no dominant features such as pyramids or round/elongated pits. Note that the irregular features shown are only three or four atomic layers deep and that the colouring exaggerates the differences in depth. Macroscopically, these surfaces have a smooth appearance and they are preferred in comparison to the surfaces showing larger shallow round pits or pyramids. The reason for this type of morphology is that the removal of HSD and VSD sites is still rather isotropic but, in comparison to the TD sites, not much more frequent. In other words, the anisotropy between pit nucleation and step propagation has been reduced, and as a consequence, no dominant features can appear, except for very small pits (holes), which are only several layers deep. The fact that the etching process becomes more isotropic with increasing coverage results from the fact that TD atoms are removed more frequently due to the increasing number of OH groups on surface, before the interaction between the OH groups at indirect second neighbours becomes relevant at higher coverages. This is the same reason why the etch rate of any orientation shows a maximum in the dependence on coverage [36].

Morphology of other orientations
In this section we consider the morphology of other surface orientations for cases A and B. It is shown that the surface features of any orientation can always be related to either exact Si(110), misaligned Si(110), exact Si(111) or misaligned Si(111), in agreement with experiment [14]. In addition to demonstrating that the most relevant morphologic features in wet etching have an atomistic origin, the most important result of this section is that the mechanism producing zigzag structures on vicinal (110) is not the same as the mechanism producing pyramidal hillocks on (100).
We start by considering the morphologies of (110) (section 4.1) and (111) (section 4.2) in the conditions that lead to the formation of pyramidal hillocks on (100) (case A). The morphologies obtained in case B are presented in section 4.3. Figure 12 shows the change in morphology as the surface orientation is gradually varied from (110) towards (100) in case A. Nosed zigzag structures [14]- [16,28] appear gradually as the misaligment is increased. At exact (110), no zigzag structures are seen at the length scale of these simulations. It may be possible that precisely aligned parallel ripples develop in real (larger) systems. For small misorientations, the zigzags are long and thin. For higher misaligments, the zigzag noses become the dominant feature as their size and density increases. The evolution of the surface is characterized by the motion of the zigzag noses from south to north along the direction parallel to the long edges. Note that if this series were continued up to a misalignment of 45 • , the morphology would be that of figure 2(f). The mechanism leading to the formation of these structures is a combination of the anisotropy of the etching process with the misorientation of the surface. The misalignment allows for the exposure of (111) surfaces, which make up the facets of the zigzag structures. The anisotropy between the removal of sites of type 3A, at the (111) facets, and of types 3B and 2B/2C, at the 110 ridges, generates the noses and their motion. Note that this mechanism differs from the stabilization mechanism that generates the pyramidal hillocks on (100), used by van Veenendaal et al [28] in order to explain also the formation of zigzags. This difference is demonstrated in section 4.3, where similar structures are shown to appear on vicinal (110) in the absence of pyramidal hillocks on (100). Figure 13 shows the changes in morphology observed as the surface orientation is varied from (110) towards (111) or (101) in case A. For any misalignment towards (111) (figures 13(a) and (b)), the morphology is characterized by rather straight steps. However, if the misorientation is towards (101), the formation of nosed zigzags becomes the dominant feature (figures 13(c) and (d)). Furthermore, as the misorientation is increased, the zigzags eventually become polygonal terraces (figures 13(e) and (f)). The reason for this different behaviour with varying orientation towards (111) or (101) can be found in the anisotropy of the etching process, combined with the misorientation of the surface. This is explained in section 4.2. Figure 14 shows the change in morphology observed as the surface orientation is varied from (111) to (100) for the conditions of case A. Note that, except for exact (111), where triangular pits are formed, the morphology of the misaligned surfaces is characterized by polygonal terraces, more closely bunched as the misorientation with respect to (111) increases. The formation of 100.23  (111) consists of (111) terraces separated by [112] steps, which are populated by stable 3B atoms, the straight steps are formed during the etching process as a result of the one-dimensional 'peeling' of the ideal-structure steps. However, for highly misoriented surfaces between (110) and (101), the ideal structure is that of (111) terraces separated by steps containing quickly removed 2B/2C sites. This leads to the appearance of the slow [112]-step family as the segments of the polygonal structures exactly in the same way as in the previous case for misoriented (111) in figure 14. Note that the morphology of the series in figures 13(a) and (b) and in figures 14(c)-(i) will be exchanged if the triangular pits on (111) are rotated through 60 • , a process that can occur with an increase in concentration or by changing the etchant [30] due to a change in the relative magnitude of the reaction rates of 3B and 2B/2C sites. This is a dramatic macroscopic effect on the surface morphology resulting from the competition between two atomistic processes. Figure 15 shows the characteristic morphologies observed under conditions producing shallow round pits on (100) (case B), instead of pyramidal hillocks (case A). In all cases (except, of course, (100)), the morphologies are similar to those of case A, and can be explained by the same mechanisms. This allows us to conclude that the mechanism of nosed zigzag formation is not controlled by the existence of stabilizing impurities, as is the case for the formation of pyramidal hillocks on (100). The mechanism leading to the formation of nosed zigzags is simply a combination of the anisotropy of the etching process and the misalignment of the orientation, as explained in section 4.1.

Etch rates
The etch rates of all considered orientations for case A follow the expected Arrhenius behaviour as a function of temperature, as shown in figure 16. As discussed in section 3.3, the etch rate of (100) experiences a large decrease at all temperatures after texturization of the surface is completed ( figure 16 (a)). Interestingly, all orientations in the (110)-(100) series, including (100) after texturization, exhibit essentially the same activation energy (≈0.43 eV) and the only notable difference is their relative vertical shifts. Note that a similar behaviour is shown in the (110)-(101) series in figure 16(b). This reflects the fact that the same atomistic processes are controlling the etching process in all these cases, as has been explicitly shown for the case of (110) and texturized (100) in [57]. Thus, such different morphologies as those in figures 12(a)-(f), 13(c)-(f) and 14(b) are all controlled essentially by the same atomistic processes. Note that the vertical shift of the lines has a geometrical origin, the dependence of the fraction of 3B atoms on the orientation of the surface [57].
In figures 16(c) and (d), a gradual increase in the activation energy occurs as the surface orientation becomes closer to (111). This is explained by the fact that, as the orientation is 100.26 changed from (100)/(110) to (111), the relative contribution of the pit-nucleation process on (111) terraces (i.e. the removal of type 3A sites) to the activation energy is increased. Figure 16(e) presents the etch rates in case A as a function of misorientation with respect to (110) or (111), depending on the case. The maximum etch rate is obtained for (100). (110) displays a sharp local minimum as a function of misorientation towards (101) and (100). However, (110) is a local maximum as a function of misorientation towards (111). (111) is a global sharp minimum. Note that the smooth curves in figure 16(e) correspond to fits of the etch rate according to the expression R ∝ sin a (α) (where α is the misorientation) for the behaviour in the regions (110)-(100), (110)-(101) and (111)-(100), and according to the expression R ∝ cos a (α) for the behaviour in region (110)-(111). Interestingly, the exponents a follow an Arrhenius behaviour with temperature, as shown in figure 16(f). This means that, by knowing the temperature dependence of the exponent, the etch rate can be interpolated for any related orientation.

Conclusions
It is shown that the most relevant surface features observed during anisotropic wet chemical etching of crystalline silicon, such as pyramidal hillocks and shallow round pits on (100), nosed zigzags on vicinal (110), polygonal and straight terraces on vicinal (111) and triangular pits on exact (111), have their origin at the atomistic scale.
It is concluded that the pyramidal hillocks on (100) are the result of local stabilization of distributed apex atoms by (metal) impurities from solution. Further work is needed in order to address the technical question of whether the impurity atoms replace the singly coordinated silicons (as assumed in this study) or the impurities simply terminate them or interact in some more complicated manner. Also, the possibility that the impurities may catalyse the re-deposition of silicon should be considered.
The temperature dependence of the density of hillocks follows an Arrhenius behaviour, in agreement with experiment. The activation energy is a complicated average of the microscopic processes that lead to singly coordinated atom disclosures and the removal of impurities, and cannot be simply associated with only one of these processes. The simulations explain the reduction in the density of hillocks with increasing concentration/coverage (in agreement with experiment) as the result of a reduction in the removal probabilities experienced by all particle types.
In the absence of the stabilizing mechanism, shallow round pits are formed on Si(100) as a result of the anisotropy between (one layer deep) pit nucleation and (isotropic) step propagation. In the simulations, the density of round pits increases with temperature, as expected, and decreases with increasing concentration/coverage, in agreement with experiments. Two alternative mechanisms for the elimination of round pits with increasing concentration/coverage are presented. In the first one, the isotropic step propagation at low coverage becomes anisotropic at high coverage, forcing the appearance of one layer deep elongated pits and preventing the formation of (deeper) round pits. In the second case, the etching process becomes more isotropic as the reaction rates of (one layer deep) pit nucleation and (isotropic) step propagation become similar for medium coverage values and, as a result, no dominant features can appear.
It is shown that nosed zigzag structures at vicinal (110) are the combined result of the misaligment and the etching anisotropy. Thus, it is concluded that the nucleating mechanisms of morphologically related structures such as pyramidal hillocks and nosed zigzags are not necessarily the same.

100.27
The simulations confirm that the formation of (one layer deep) triangular pits on exact Si(111) and of polygonal (saw-shaped) and straight steps between the terraces in vicinal Si(111) is controlled by the relative rate of [121] and [121] step-propagation and depends on the misorientation of the surface with respect to Si(111).