Orientation-dependent etching of silicon by fluorine molecules: A quantum chemistry computational study

Orientation-dependent etching of silicon has been used for texturing of silicon surfaces via dry and wet etching processes to produce black silicon.^1,2 The etch rate in the (111) direction is much slower than in the (100) and (110) directions, and the resulting anisotropic etching process can be used to create nanostructures. Plasmaless atmospheric dry etching (ADE) commercialized by Nines Photovoltaics is a cost-effective method of nanoscale silicon surface texturing. A schematic of this process is shown in Fig. 1, illustrating the experimental etching results from Ref. 3 at 3.33% F₂ in air and substrate temperatures of 200 and 300 °C. In the ADE process implemented by Nines Photovoltaics, a silicon wafer is transported on a heated conveyer through a reaction chamber containing the F₂ gas. The substrate temperature, F₂ concentration, and duration of exposure to F₂ all control the etch process. Kafle et al.^3–6 showed that at temperatures higher than 200 °C, nanoscale pits are formed on a partly oxidized Si(100) surface during plasmaless etching. Nucleation of the pits due to anisotropic etching begins at reactive sites where the oxide layer is absent. The slope of the pit walls is defined by the angle between the Si(111) surface and Si(100) surface; these planes intersect at an angle of about 55°, as shown in Figs. 1(b) and 1(d). The cone-shaped nanostructures shown in Fig. 1(b) gradually disappear with increasing temperature [see Fig. 1(c)] as the etching process becomes more isotropic. It is schematically shown in Fig. 1(d) how the anisotropic etching profile appears on the Si(100) surface due to the slow etching of Si(111). Our focus is limited to the analysis of anisotropic etching on a molecular level, using DFT (density functional theory) simulations of reaction pathways and investigating the suitability of a reactive bond order (REBO) molecular dynamics potential for simulating this anisotropic etching process. We did not study the effects of F₂ concentration or process duration on etch anisotropy during ADE processing; nor have we examined the process of pit nucleation on the partly oxidized surface. It should be noted that a similar effect, where silicon is slowly etched in the (111) direction compared to the (100) direction by Cl₂/Ar inductively coupled plasma, was observed by Du et al.⁷ 

To identify the reasons for anisotropy, we performed quantum chemistry modeling of F₂ dissociative chemisorption on silicon surfaces. In our modeling, we assume that the etching of silicon by F₂ proceeds on the F-terminated surfaces and that the rate of silicon etching is determined by the rate of F₂ chemisorption. The calculated reaction rates of F₂ dissociative chemisorption on Si(100) and Si(110) match the rate of silicon etching by F₂ gas as measured by Mucha et al.⁸ (see Sec. II D). In addition, our calculations reproduce the experimental trend that F₂ gas etches the Si(111) surface more slowly compared to other surface orientations as illustrated in Fig. 1(d). The charge density analysis (described in Sec. II C) showed that passivating F atoms attract electron density from the uppermost Si atoms and that the amount of positive charge on an individual Si atom increases with the number of bonded passivating F atoms. Si(100) and Si(110) surfaces are passivated by a larger amount of F atoms per surface Si atom compared to Si(111) (Sec. II A). F₂ molecules approaching the surface become negatively charged and, as a result, are more strongly attracted to Si(100) and Si(110) surfaces than the Si(111) surface. The analysis shows that the etching limiting step is F₂ dissociative chemisorption on silicon, and their activation barriers (E_a) for Si(100) and Si(110) are approximately equal, E_a [Si(100)] ≈ E_a [Si(110)], and are small compared to E_a [Si(111)] (Sec. IV).

Additionally, we performed the classical molecular dynamics simulation of etching. It was determined that the REBO potential that was typically used for this chemistry reasonably represents chemisorption on Si(111) but did not accurately represent the reaction path for F₂ dissociative chemisorption on Si(100) and Si(110) surfaces, overestimating the required activation energy. As a result, the anisotropic orientation-dependent etching process that creates conelike nanostructures cannot be simulated using this potential without further reparameterization. The results can be found in Sec. III.

As is well known, bulk silicon crystallizes into a diamond cubic crystal structure in which each atom has four nearest neighbors forming the sp 3 hybridized tetrahedral structure. The structures of fluorinated Si(100), Si(110), and Si(111) surfaces and the top 6 silicon atoms, so-called Si₆ chairs, on these surfaces are shown in Fig. 2. The reaction path of F₂ dissociative chemisorption passes via the transition state (TS), where an F₂ molecule adsorbed on the Si₆ chair forms a reaction center with attacked Si–Si. In other words, only two atoms in the Si₆ chair are partaking in the reaction, and the stoichiometry of the reaction centers is F₂..Si₂F_n. The reaction centers have different “n” depending on the surface orientation. Namely, there is an F₂..Si₂F₂ reaction center on Si(100) and Si(110) surfaces and an F₂..Si₂F₁ reaction center on the Si(111) surface as illustrated in Fig. 2. Note that the reaction centers shown in Fig. 3 are coupled with possible reaction channels A and B: one F atom of the F₂ molecule bounces off the surface yielding a free F atom or both F atoms are adsorbed without the formation of free atoms. Both channels are illustrated in Fig. 3 for all Si₆ chair configurations.

It should be noted that the Si(100) surface can undergo reconstruction to minimize the energy of the top layer.⁹ Therefore, in Sec. II C, we also consider F₂ dissociative chemisorption on F-terminated Si(100) reconstructed surfaces as well. We showed that the reconstructed F-terminated Si(100) surface is more reactive (less stable) with F₂ molecules than the unreconstructed Si(100) surface. Thus, we assumed that the etching of silicon in the (100) direction is not accompanied by reconstruction and proceeds on the unreconstructed Si(100) surface.

As mentioned in Sect. II A, the chemisorption of F₂ proceeds via two competitive reaction channels known as chemisorption with single atom abstraction (channel A) and with two atom adsorption (channel B). Channel A requires four elementary steps to etch one silicon atom, while channel B requires only two, as shown in Fig. 4. The first step is the same for both reaction channels A and B. The divergence between channels A and B appears when the TS goes down to the reagent valley, resulting in two different products (see Sec. II C). The importance of channel A is that free F atoms are released, which should affect the etch profile. These highly reactive free F atoms can also etch the silicon surface, but rather than the anisotropic, orientation-dependent etch that is possible with molecular fluorine, atomic F etches Si(100), Si(110), and Si(111) with the same rate. Thus, chemisorption that proceeds via channel A causes the etch profile to become more isotropic.

Thus, we performed quantum chemistry modeling and found TSs of the first step in F₂ chemisorption on Si(100), Si(110), and Si(111) surfaces. Other steps and evolution toward products in channels A and B are not considered. However, the calculated reaction barriers of step 1 match the measured activation energies for silicon etching by molecular fluorine (see Sec. IV). This agreement confirms our assumption that step 1 is the rate-determining step of the etching process.

Quantum chemistry modeling was performed using broken symmetry orbitals calculated by the unrestricted u-B3LYP DFT functional in the gaussian 16¹⁰ software package. The Lanl2dz basis set and pseudopotential for Si atoms were used in our simulations. All structures were visualized by Chemcraft software.¹¹ 

An Si₃₂F₃₂ cluster was used to mimic the structure of fluorinated Si surfaces. The free bonds of external Si atoms were closed by fluorine atoms. The reaction center geometries of four transition states (TSs) corresponding to the reactions of F₂ dissociative chemisorption on an Si(100) surface with and without surface reconstruction, on an Si(110) surface, and on an Si(111) surface are shown in Fig. 5. In all transition states considered here, the F–F bond is perpendicular to the Si–Si bond. On the reaction path after the transition state, an intermediate biradical structure is observed where singlet and triplet states are overlapping, as shown in Fig. 6. Natural population analysis was performed to calculate atomic charges. The charges on the Si₂F_n reaction centers on the Si₃₂F₃₂ cluster before the reaction and in TS can be found in Fig. 7. (The Si₆ chair of the Si₃₂F₃₂ cluster reacting with an F₂ molecule is shown in the displayed structure, while other atoms were omitted to simplify the picture.) The charge distribution is quite complex; the total charge of the reaction center is not zero, but the total charge of the whole F₂..Si₃₂F₃₂ molecular system is zero. As shown in Fig. 7, the charge on the Si atoms of the F₂..Si₂F_n reaction centers increases slightly as F₂ molecules approach the Si–Si bond due to attraction of electron density by F₂ molecules. However, the majority of positive charge on these Si atoms is due to electron displacement to F substitutes that form a part of the reaction center not due to electron displacement to F₂ molecules. Thus, the reactivity of the surface is determined by the surface stoichiometry, which, in turn, determines the surface charges and the reaction barrier for F₂ dissociation.

The presence of the biradical structure in Fig. 6 is consistent with dissociative chemisorption via single atom abstraction as proposed by Tate et al.,^13,14 in which only one atom of F₂ forms an Si–F bond, while the other F atom desorbs into the gas phase [channel A in Fig. 4 and product A in Fig. 6]. Channel B in Fig. 4 and product B in Fig. 6 correspond to the reaction channel of dissociative chemisorption, where both atoms of the F₂ molecule adsorb on the silicon surface resulting in two Si–F bonds. Note that molecular dynamic simulations by Carter et al.¹⁵ showed that both reaction channels yielding product A and product B are present, and they are the dominant reaction channels on a clean, unpassivated Si(100) surface. We expect that the reaction channel with single atom abstraction decreases the selectivity of orientation-dependent etching because one of the F atoms of the F₂ molecule (which initially desorbs into the gas phase) is redeposited on the substrate and nonselectively etches Si(100), Si(110), and Si(111) surfaces.

The results are presented in Fig. 8, which shows that our calculations give good agreement with experimental data if we assume that the exposed surface in the experiment is the Si(100) facet or the Si(110) facet. [We arbitrarily chose the density of n(F₂) to be equivalent to 1 Torr partial pressure in order to present the result in convenient units, monolayers per second.] Note that the effect of atomic F released by the single atom abstraction mechanism is not considered in our calculation of the etch rate in Eq. (7) and Fig. 8.

Pullman et al.¹⁶ experimentally studied dissociative chemisorption of F₂ on an Si(100)(2 × 1) surface saturated with a single monolayer of fluorine by exposing the fluorine-saturated surface to supersonic F₂ beams of variable energy. They found that no reaction occurs below 0.16 eV of incident energy. Note that this experiment aimed to study the reaction of F₂ dissociative chemisorption without subsequent etching. Therefore, we assume that the studied reaction proceeded on initially relaxed Si(100) with surface reconstruction. Indeed, our calculated reaction barrier for F₂ dissociative chemisorption on the reconstructed Si(100) surface (0.13 eV) is close to the measured threshold value of F₂ incident energy (0.16 eV). In other words, the measured threshold value of F₂ incident energy is determined by the barrier of F₂ dissociative chemisorption on the reconstructed surfaces.

The reaction of molecular fluorine with an Si(111) surface was studied by Tatsumi and Hiroi¹⁷ in an ultrahigh vacuum reactor. They found that the effective activation energy for the surface reactions of molecular fluorine leading to Si(111) surface etching at temperatures below 580 °C was 0.61 eV. According to our calculations, the reaction barrier for F₂ dissociative chemisorption on Si(111) is 0.57 eV.

We also performed classical molecular dynamics (MD) simulations to deepen our understanding of crystallographic orientation-dependent etching of Si by F₂. In these simulations, Si–F interactions are described using a reactive empirical bond order (REBO) potential. These types of force fields are commonly used to simulate covalently bonded materials (such as Si) as well as their reactions with a variety of gaseous species. We use the REBO potential parameterized by Humbird and Graves.¹⁸ 

In order to utilize the MD code in an efficient manner, we simulate the etching process by a series of “impact simulations.” These simulations consist of a semi-infinite Si slab with a vacuum space above the slab in the z direction. The Si slab consists of Si atoms in a diamond cubic lattice. Periodic boundary conditions are imposed in x and y directions to make the slab semi-infinite. The bottom two layers are fixed to prevent drift of the slab due to incoming F₂ impacts. During an “impact simulation,” an F₂ molecule is placed in a random position in the vacuum space above the Si slab. The height of the molecule is chosen such that it is just outside the interaction distance of the Si slab. The velocity components of F₂ are randomly sampled from the Maxwell–Boltzmann distribution at 300 K. We only consider molecules with a negative z component of velocity, so the molecule will impact the Si slab. Next, a microcanonical ensemble (constant number of atoms, constant volume, and constant total energy) MD simulation is run lasting approximately 2 ps. After the trajectory is complete, species that have been etched or sputtered are deleted from simulation. In addition, species that are weakly bound to the surface are also deleted. After this, a Berendsen thermostat¹⁹ is applied to the simulation cell to remove excess energy introduced by the incoming species and to maintain the temperature of the Si slab. More information on the thermostat procedure and product deletion routines can be found in previous studies.^20,21 This procedure is repeated until the desired fluence of F₂ is reached. We consider Si slabs with different crystal facets facing the vacuum space to study the different etching behavior of the facets. The surface area for the Si slab with the (100) surface exposed (both unreconstructed and reconstructed) is about 472.5 Ų. For the slab with the (110) surface exposed, the surface area is about 626.4 Ų. The surface area for the Si slab with the (111) surface exposed is roughly 613.8 Ų. The depth of each cell is large enough such that no F₂ molecules will interact with the fixed layer. The lattice constant in all cells corresponds to diamond cubic Si at 300 K. Even though we consider temperatures higher than 300 K, the lattice constant remains the same in our simulations. In reality, the crystal will expand as the temperature is raised. However, because the thermal expansion of Si is very small (the lattice constant increases by about 0.5% from 298.5 to 1513.2 K²²), we ignore this expansion in our simulations.

Our classical MD simulations did not show silicon etching at temperatures below 900 K. (F₂ bounced off of the silicon surfaces every time.) In order to explain the absence of etching, we performed intrinsic reaction coordinate (IRC) DFT calculations of F₂ dissociative chemisorption on Si₃₂ clusters for all TSs shown in Fig. 3. IRC calculations follow the steepest descent path from the TS to the reagent and product valleys. The same reaction pathways were recalculated using the REBO classical MD potential. As shown in Figs. 9(a)–9(c), very high barriers appear on the (100) and (110) reaction paths recalculated by the classical MD potential, which do not appear on the DFT-calculated reaction paths. The (111) reaction path, on the other hand, is quite accurately reproduced by the REBO potential.

The lack of etching observed in molecular dynamics simulations at temperatures below 900 K can be explained by high reaction barriers. The reactions are very slow below 900 K and, as a result, cannot be simulated by classical MD in a reasonable amount of time. In the case of Si(100) and Si(110) surface etching, the barriers are overestimated and additional fitting of the potential would be required for Si + F₂ chemistry in order to correct this. Without reparameterization, the MD potential would be unable to simulate the anisotropic etching described earlier in the paper. The reaction barrier for F₂ dissociative chemisorption on Si(111) calculated by the REBO potential is very close to the DFT value (0.61 vs. 0.57 eV). However, the inverse reaction rate (timescale of the reaction) is still higher than 1 ns at temperatures below 900 K [see Fig. 6(d)]. The etching probabilities calculated by classical MD and the probabilities of F₂ dissociative chemisorption calculated by transition state theory (TST) can be found in Fig. 10.

As our DFT calculations show, the dissociative chemisorption of F₂ on Si(111) proceeds over a significantly higher barrier than the dissociative chemisorption on Si(100) and Si(110) surfaces. This results in a lower reaction probability [Fig. 10] and, hence, a lower etch rate for the Si(111) surface. We also demonstrate that etching of the (100) surface with reconstruction is faster than etching on the (100) surface without reconstruction. We assume that the removal of surface layers by fluorine proceeds without reconstruction. Thus, it is equally probable that the F₂ molecule dissociates on unreconstructed (100) and (110) surfaces, so etching in the (100) and (110) directions proceeds at the same rate. The low etch rate of the Si(111) surface, by comparison, explains the orientation-dependent etching that was used by Kafle et al.^3–6 for texturing silicon surfaces during black silicon production.

We attribute the high etch rate of Si(100) and Si(110) to the fact that these surfaces are passivated by a larger number of F atoms per surface Si atom compared to Si(111) surfaces. The electron density is displaced from Si atoms to more electronegative passivating F atoms, resulting in a partial positive charge on Si atoms [Fig. 11]. (This effect is known as the inductive effect in organic chemistry). The amount of positive charge on a surface increases with “n,” and as a result, Si(100) and Si(110) surfaces accumulate larger amounts of positive charge relative to the Si(111) surface. Thus, the reaction barrier for dissociative chemisorption of molecular fluorine is smaller for Si(100) and Si(110) surfaces compared to Si(111) surfaces [Fig. 11(b)]. An F₂ molecule approaching the silicon surface becomes negatively charged, attracting electron density from surface Si atoms. As a result, F₂^δ− is more strongly attracted to Si(100) and Si(110) surfaces than to the Si(111). This leads to observed fast etching in (100) and (110) directions and slow etching in the (111) direction.

Our simulations show that the activation barriers of F₂ dissociative chemisorption are 0.31 eV on Si(100), 0.35 eV on Si(110), and 0.57 eV on Si(111). Finally, we collect calculated parameters for rate constants and probabilities of F₂ dissociative chemisorption on silicon surfaces in Table I.

Orientation-dependent etching of silicon via molecular F₂ has been studied using quantum chemistry methods. Activation barriers, reaction rates, and probabilities for F₂ dissociative chemisorption on (100), (110), and (111) surfaces were calculated. It was shown that the reaction barrier of F₂ dissociative chemisorption is a function of the surface stoichiometry, which, in turn, determines the amount of positive charge on surface Si atoms. As a result, the etching of Si(111) is significantly slower than (100) and (110) surfaces. Our simulation results are consistent with previously published experimental data and explain the mechanism of orientation-dependent etching of Si by F₂ gas.

Additionally, we describe two reaction channels of F₂ dissociative chemisorption. We expect that the desorption of F atoms during F₂ dissociation is significant during the etching process. This is consistent with dissociative chemisorption via single atom abstraction, which was proposed previously. We expect that the reaction channel via single atom abstraction causes nonorientation-dependent etching with the isotropic profile.

We also used molecular dynamics simulations to investigate the suitability of a REBO potential for modeling silicon etching via molecular F₂. While the REBO potential reasonably described F₂ dissociative chemisorption on Si(111), it predicts unrealistically high reaction barriers for dissociative chemisorption on Si(110) and Si(100). This results in very slow etching rates and an inability to reproduce the anisotropic etching mechanism studied in this work. In the future, molecular dynamics modeling of this process will require reparametrized potentials that accurately describe orientation-dependent reaction barriers.

The research described in this paper was conducted under the Laboratory Directed Research and Development (LDRD) Program at Princeton Plasma Physics Laboratory, a national laboratory operated by Princeton University for the U.S. Department of Energy under Prime Contract No. DE-AC02-09CH11466. The United States Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

In addition, this research used computing resources on the Princeton University Adroit Cluster and Stellar Cluster. The authors gratefully acknowledge discussions about the classical MD code with David Humbird (DWH Process Consulting).

The authors have no conflicts to disclose.

Omesh Dhar Dwivedi: Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Yuri Barsukov: Conceptualization (equal); Investigation (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Sierra Jubin: Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Joseph R. Vella: Investigation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Igor Kaganovich: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and from the corresponding author upon reasonable request.