Revealing the Chemical Reaction Properties of a SiHCl3 Pyrolysis System by the ReaxFF Molecular Dynamics Method.

The pyrolysis kinetics of SiHCl3 and its reaction mechanism are essential for the chemical vapor deposition process in polysilicon industries. However, due to the high temperature and lack of in situ experimental detection technology, it is difficult to carry out experimental research on the pyrolysis kinetics of SiHCl3. In this work, reactive force field molecular dynamics simulations of SiHCl3 pyrolysis were performed to investigate the effect of temperature on the pyrolysis kinetics of SiHCl3 at the atomistic scale in a wide temperature range (1000-2000 K). The lumped Si clusters containing more than five Si atoms tended to appear at the later period of the reaction under a temperature lower than 1300 K, some of which even possessed polycyclic structures; nevertheless, small ones with less than two Si atoms such as SiHCl2 and HCl tended to emerge under a high temperature. The changes of partial energy terms with time evolution under various temperatures were proved to be rooted in the distribution of intermediates based on the momentary simulation period. In general, the reaction network at a low temperature was more complicated than that at a high temperature, resulting from the fact that more chemical events and intermediates came into existence, and the maximum number of Si atoms in one single molecule/radical was observed under a low temperature than that under a high temperature. As to the variation of SiHCl3 with the progress of the reaction, the linear fitting tendency disappeared under the temperature above 1300 K, which changed in fluctuation with the further elevation of temperature, elucidating the fact that SiHCl3 can act as a product and not just as a reactant to participate in elementary chemical events frequently.

Introduction

A great deal of efforts have been made over the past decade in facing the global fossil fuel crisis and environmental problems, leading to a rapid increase in demand for renewable energy, especially solar energy. As a consequence, investigations relevant to high-purity polysilicon production with high efficiency and low cost are of great significance due to the important role it plays in accelerating polysilicon production since high-purity polysilicon is an indispensable raw material for the preparation of photovoltaic solar cells.[1] In general, high-purity polysilicon is mainly produced by the thermal decomposition and reduction of trichlorosilane (TCS, SiHCl3) inside a H2-rich chemical vapor deposition (CVD) reactor with a temperature over 1000 °C,[2,3] accompanied with the formation of various byproducts such as tetrachloride (STC, SiCl4), dichlorosilane (DCS, SiH2Cl2), and so on, making it too complex to investigate the possible chemical events and intermediates involved in the CVD process.

Many efforts have been made in calculating the TCS reaction system in theory effected by the motivation from industry upgrading, such as the first-principles calculation, whose advantages lie in analyzing complex reaction systems, intermediates of fast reactions, and elementary reactions with low proportions. For example, the hydrochlorination reaction with the assistance of Cu–Si alloy catalysts during polysilicon production was investigated by the plane-wave-based density functional theory (DFT) method,[4] in which the adsorption behaviors of STC, TCS, DCS, SiCl2, HCl, and H2 were compared at two types of catalyst surfaces. Dissociative chemisorption reactions were more likely to happen at the Si-terminated surfaces for the exposed Si atoms on those rough surfaces, while no molecular adsorption occurred on the Cu-terminated surfaces owing to a flat geometry and a uniformly distributed charge. Similarly, the dissociative chemisorption of H2, HCl, SiCl2, DCS, TCS, and SiCl4 on the periodic slabs of the Si(100)-c(4 × 2) surface was studied using the DFT method,[5] elucidating that the strongest binding energy was observed in HCl and SiCl2 among all these molecules. Aside from the chemical adsorption process, pervasive use of chlorosilanes and their various derivatives has led to extensive research associated with their reactions between free radical specials. For instance, the ab initio direct dynamics method was used to explore the reaction path dynamics and theoretical rate constants for the hydrogen abstraction reaction of SiH3Cl + H* → SiH2Cl + H2 as well as the influence of chlorine substitution on the reaction reactivity trends with the discussion of many related factors, including the variation of bond lengths, vibrational frequencies, potential energies, and the total curvature along the reaction path.[6] As a result, the three-parameter fit for the reaction rate constants was acquired as k = 1.5 × 106T2.40 exp(−1134/T) cm3 mol–1 s–1 over a wide temperature range from 200 to 3000 K. However, these studies only allow the simulation of a specific reaction or a class of reactions but did not provide an atomistic description of the extremely complex and complete reaction process involved in many chemical events.

Besides, the comprehensive reaction network in the gas phase composed by a large number of elementary reactions was extensively researched. Based on the quantum chemical calculations and transition state theory, the reaction mechanisms of SiH4–Cl (n = 0–4) in the gas phase were utilized in the system containing 20 species and 117 elementary reactions, whose rate coefficients were acquired through calculating the thermodynamic properties.[7] Similarly, the dominant chemistry active in the gas phase for the CVD process of Si from TCS was explored systematically to compare two different routes for the gas-phase reactivity, namely, disilane mechanism and radical mechanism, with the creation of a reaction model containing 26 reactions and 16 chemical species.[8]

Although the reaction models including single or few chemical events being used in the above valuable research studies allow for the simulation of a few reactions in a complicated reaction pathway, it was difficult to form an integrated description of all the possible chemical events that occurred in one process since these chemical events were pre-defined and number-limited. At this point, a detailed mechanism and product distribution relevant to all the possible chemical events in the CVD process of Si from TCS were still absent. In general, the extremely complex CVD process of Si from TCS was constituted by the H2 environment, gas-phase reaction, and the interface reaction, both of which were significant and indispensable to get a good knowledge of the reaction mechanism of the CVD process in the micro-scale as a whole. In other words, the complicated process could often be decomposed into multiple relatively simple processes. As one of the vital sections and the most basic step of the CVD process, the investigation of the pyrolysis process of TCS could provide help to precisely distinguish the role of a H2 environment in the gas-phase reaction on one hand and quantitatively estimate the influence of the interaction between the Si-base and TCS on the generation of monatomic silicon in the surface reaction on the other hand. Besides, the pyrolysis kinetics of SiHCl3 acquired in the present study could be transported to other investigations and provide great convenience for these investigations where H2 was not allowed for some reasons. Additional simulations could be performed in our future work to study the complete gas-phase reaction with TCS and the H2 environment and the surface reaction in the gas phase and the solid phase with TCS, H2, as well as the Si-base in the micro-scale. Thereafter, the TCS pyrolysis system in the gas phase was investigated in the present work by employing the reactive force field molecular dynamics (ReaxFF MD) method in order to acquire the dynamic properties and product distribution of the TCS pyrolysis process. It had been demonstrated that a quantum mechanical (QM)-derived ReaxFF MD method could be applied to explore complicated chemical processes over long time scales for large systems such as elucidating the chemical events relevant to the explosion of high-energy materials,[9] pyrolysis and combustion,[10−13] and interactions between crystals and crystal surfaces.[14,15] The sole derivation of the ReaxFF force field parameters from QM allowed for the direct application of the ReaxFF MD method to novel systems that might not have been or could not be extensively studied experimentally.[16]

Specifically, a series of simulations associated with the TCS pyrolysis process over a wide temperature range between 1000 and 2000 K were presented to understand the pyrolysis kinetics of TCS at the atomistic scale. Not only were all the possible intermediates presented in detail, but also the close relationship of the dynamic evolution between intermediates and energy terms was explored and evaluated in this work, all of which were very helpful for understanding the CVD process of Si from TCS.

Results and Discussion

Analysis of the Main Intermediates’ Distribution

In summary, the TCS pyrolysis system was simulated by the ReaxFF MD method to obtain a good knowledge of the detailed pyrolysis mechanisms, including the dynamic and fundamental products’ distribution and the pyrolysis temperature dependent on the dynamic evolution of the intermediates.

At first, part of the products’ distribution and reaction channels were reported and discussed, aiming at validating the accuracy of the simulation results. One kinetic mechanism constituted by 26 reversible reactions was provided in the research that was devoted to distinguishing the dominant gas-phase chemistry active during the CVD process of Si from TCS through computing the rate constant of each reaction, which consist of the following 16 chemical species: H2, HCl, SiHCl3, SiCl4, SiH2Cl2, SiCl2, Si2Cl6, Si2HCl5, HCl2SiSiCl2H, H2ClSiSiCl3, Si2Cl4, SiCl3, SiHCl, SiHCl2, H, and Cl.[8] In the present work, intermediates of HCl, SiHCl3, SiCl4, SiCl2, Si2HCl5, HCl2SiSiCl2H, SiCl3, SiHCl2, H, and Cl were produced in the system under 1800 K, repeating 10 of these 16 chemical species. In contrast, there were 10 products that were identical to these 16 chemical species when the temperature reduced to 1600 K, which just replaced SiCl2 by SiH2Cl2 in comparison with that under 1800 K. Similarly, the number of intermediates that duplicated the above 16 chemical species was 11 with the continuous declaration of temperature until the lowest temperature of 1000 K, namely, the newly increased Si2Cl6 compared to that under 1800 K. Indeed, our simulations under diverse temperatures just reproduced approximately 70% of these 16 chemical species but should be considered as good agreement given the uncertainties of the calculated data. Nevertheless, these 16 chemical species did not cover all of the possible intermediates. On the other hand, the recurrence of some products with the distinctive structures in the present work could also play a role in proving the present simulation results to some extent, such as the transition states located for the reaction: Si2HCl5 → SiCl2 + SiHCl3, namely, , which was captured in the system under 1900 K. Besides, the present simulations were able to regenerate the dominant reaction channels provided in the other calculation work, suggesting the credibility of the present work. For instance, the reaction: SiHCl3 → SiCl2 + HCl was regarded as the significant one for initiating the decomposition of SiHCl3 in the investigation of the gas reactivity of chlorosilanes[8] and the SiH4–Cl (n = 0–4) reaction mechanisms for the polysilicon production process.[7] In fact, the present simulations partitioned this step into two reverse steps: 2SiHCl3 ⇌ Si2HCl5 + HCl, followed by Si2HCl5 ⇌ SiCl2 + SiHCl3. In conclusion, the precision of the present work could be verified by all of these comparisons.

In order to confirm the reasonableness of the 1000 ps simulation time employed to obtain the pyrolysis kinetics of TCS, the number of reactants and intermediates with the time and temperature is shown in Figure , where the intermediates containing n silicon atoms were lumped as Si species. As displayed in Figure a, the number of TCS molecules decreased with time at 0–600 ps and at 1000 and 1300 K and fluctuated with time in the following process. When the temperature increased to 1600 and 2000 K, the TCS molecules experienced slight decreasing trends at 0–200 ps, and their numbers were almost unchanged at 200–1000 ps. The number evolution of TCS molecules suggested that the 1000 ps simulation time was long enough to reach the equilibrium of TCS consumption and generation in the simulated temperature range.

Figure 1

Number variation of species with the evolution of time and temperature obtained in 1000 ps simulation of SiHCl3 pyrolysis system: (a) SiHCl3 molecules, (b) Si2 species, (c) Si3 species, (d) Si4 species, (e) Si0 species.

Temperature was an important factor to influence the process. As shown in Figure b–d, the production of the species in Figure at 2000 K was lower than that below 2000 K, illustrating the adverse effect of extremely high temperature on the generation of specific species. In addition, the number variations in the trend of Si2 species > Si3 species > Si0 species > Si4 species elaborated the easier creation of the lumped Si2 species than that of Si3 and Si4 species. There was essentially no observation of Si4 species during the dynamic evolution of the pyrolysis system under 2000 K. Consequently, it was more likely that the difficulties for the creation of relatively large molecules (i.e., the lumped Si3 and Si4 species) from the polymerization of small ones gradually grew with the increase of the pyrolysis temperature, relevant to the unsuccessful collisions between small molecules. From the number variation of intermediates observed during the simulation, the lumped Si2 species were expected to be the primary intermediates compared with other intermediates and would play a crucial role in the micro-reaction mechanisms of TCS pyrolysis. Additionally, the lumped Si3 and Si4 species were obtained in the late period of the pyrolysis process (Figure c,d).

In accordance with the lumped Si species discussed above, the Si0 species was due to the class of intermediates without Si atoms, including hydrogen radicals (H*), chlorine radicals (Cl*), HCl molecules, and so on. However, Cl2 molecules were basically not generated during these simulations. One of the possible reasons that existed was the low concentration of Cl* to make a successful encounter between two Cl* to further form a Cl2 molecule. The other possible reason was that the parameters of the ReaxFF applied in the present study had not been optimized for the generation of the Cl2 molecule. Since the generation of the Cl2 molecule was not the key research content of this study, the absence of the Cl2 molecule did not play a part in the further analysis of the simulation results. Some relevant trials would be done in our future work to focus on the possibilities and rationality for the production of the Cl2 molecule. As shown in Figure e, the generation of the lumped Si0 species preferred to happen at a relatively later period of the pyrolysis process with the benefits from the low temperature. In addition, for each simulation system, the detailed intermediates embodied in the Si0 species were not identical but differed with the pyrolysis temperature. Specifically, the Si0 species comprised H*, Cl*, and HCl molecules in the system under 1300, 1600, and 2000 K, while more kinds of H*, Cl*, HCl, HCl2, and H2Cl2 intermediates constituted the lumped Si0 species being generated in the system under 1000 K. The inference that the reactive environment was more complex at a low temperature than that at a high temperature could be drawn owing to the fact that much smaller molecules and simple radicals without a silicon atom were produced in the pyrolysis system under a low temperature, which were collectively referred to as the Si0 species. Even though all possible reactions were available with the usage of the ReaxFF MD method, only the reaction pathways that were most kinetically and thermodynamically plausive would be ensued during the simulations, resulting in the dependence of the chemical events on the pyrolysis temperature. In conclusion, there was a significant influence of the pyrolysis temperature not only on the type but also on the number of intermediates created during the simulation, especially for the lumped Si0, Si2, Si3, and Si4 species.

Apart from the lumped Si4 species, the phenomenon of silicon clustering was even observed during the simulations, playing a significant role in investigating the properties of the TCS pyrolysis process. The evolution curve in Figure shows that the change of the total number of chemical reactions and total number of intermediates as a function of temperature were in agreement with each other, all of which basically build up over the reduction of the pyrolysis temperature. However, as to the total number of intermediates, the degree of pyrolysis temperature dependence became weaker at a high temperature above 1300 K. For example, the total number of intermediates obtained at 1300 K only accounted for 10% of that acquired at 1000 K, resulting from the small temperature interval of 300 K. It seemed that more kinds of molecules and radicals were generated by the enhancement of the total number of chemical reactions during the simulation. In addition, diverse silicon clusters emerged at all temperatures accompanied by the maximum number of silicon atoms in a single molecule ranging from 3 to 17. Specifically, the largest silicon aggregates gently changed with reducing temperature, but two abrupt jumps in the size of the silicon cluster appeared at low temperatures between 1300 and 1000 K. At this point, two distinct mechanisms were illustrated by our simulations, namely, one dominated by large silicon clusters beneficial for low temperatures and a second in the generation of small molecules by most Si atoms in favor of high temperatures, which were in agreement with the occurrence of C aggregates in the presence of the thermal decomposition of RDX [cyclic-[CH2N(NO2)]3].[9] To entirely explore the relatively slow course of silicon aggregation, the simulations in larger systems and longer timescales might be indispensable.

Figure 2

Total number of chemical reactions and intermediates and the maximum number of Si atoms in a single molecule as a function of temperature in an ReaxFF MD simulation of the SiHCl3 pyrolysis process.

To explore the relationship between intermediates and reaction pathways, some typical 2D structures of the molecules containing the maximum silicon atoms, which were produced in the system under 1000, 1300, 1600, and 1900 K, are described in Figure . Similarly, the composition of these 2D structures varied with temperature, and a low temperature was extremely conducive to building up the complexity of these 2D structures, especially for the creation of a Si17H12Cl45 molecule in the system under 1000 K. The polycyclic structures in the Si17H12Cl45 molecule (1000 K) could be attributed to the emergence of abundant Si–Cl–Si bonds, functioning as one bridge to join relatively small molecules into large even enormous ones. It might be rooted in the fact that the energies of 87.4 and 98.5 kcal/mol were necessary for the scission of Si–H and Si–Cl bonds, respectively, suggesting that the Si–H bond was easier to break than the Si–Cl bond, further illustrating the high reactivity of the H atom than that of the Cl atom.[7] Therefore, it was the stable Si–Cl–Si bonds and not the Si–H–Si bonds that enhanced the connectivity of Si clusters. Although the sophistication of the molecules shown in Figure could not be arranged in the same level for comparison, the bridge of the Si–Cl–Si bonds appeared in the molecules of Si4H3Cl11 (1300 K), Si5H5Cl14 (1600 K), and Si3H3Cl9 (1900 K), irrespective of the distinct molecular structures such as the linear form and the polycyclic modality. In some sense, the restriction of the thermal rate of atom motion rooting in the low energy of the system under a low temperature in comparison with a high temperature directly leads to the relatively small distance between atom pairs and enhancement of the intermolecular attractive force at the same time, favoring the formation of new bonds. In addition, compared with the Si–H atom pair, it was the Si–Cl atom pair that prior to form bond under the same condition resulting from its smaller bond order (BO) cutoff (0.55) than that of Si–H atom pair (0.6). Even though the BO cutoff of Si–Si (0.3) was much smaller than that of Si–Cl, the difficulty of the bond compression of Si–Si suffered from its less population and larger space resistance. The above discussion provided an explanation for the emergence of the extremely complex structure and the large population of Si–Cl–Si bonds, as shown in Figure . Once again, these phenomena indicated the complexity of the reaction environment at a low temperature than that at a high temperature. Accordingly, the augmenting of the pyrolysis temperature was responsible for weakening the possibilities of the occurrence of side chemical reactions and the generation of byproducts just when the temperature varied in the range below 1300 K, while the further increase of the pyrolysis temperature could not help to reduce these potentials to a great extent if the temperature was higher than 1300 K.

Figure 3

One typical 2D structure of the largest intermediate at different pyrolysis temperatures: (a) 1000, (b) 1300, (c) 1600, (d) 1900 K.

Furthermore, the sensitivity of the consumption of the reactant TCS to the pyrolysis temperature was observed during the simulations. Figure shows the time evolution of the population of TCS molecules within the timescale of 1000 ps for isothermal simulation, where the fitting of one straight line was presented under 1100 K (Figure a), while the non-monotonous and alternating variation tendency between increase and decrease appeared under 1900 K (Figure b). In fact, the transition point of this linear fitting trend into a non-monotonous one took place at the temperature of 1300 K. In other words, the variation of the population of TCS molecules with time evolution grew similar to the curve shown in Figure b at a high temperature above 1300 K. Additionally, the significant higher consumption of TCS molecules was obvious under a low temperature in comparison to a high temperature. These differences occurred since the reaction pathway of the clusterization of small molecules into large molecules by the formation of more stable Si–Cl–Si bonds was thermodynamically beneficial at a low temperature rather than at a high temperature.

Figure 4

Dynamic evolution of the number of SiHCl3 molecules at different pyrolysis temperatures: (a) 1100, (b) 1900 K.

It was expected that the type of chemical reactions involved in TCS molecules was closely relevant to the pyrolysis temperature. To get a better knowledge of this point of view, the changes of WTCS with the various temperatures studied are shown in Figure , which were calculated with respect to the proportion of elementary reactions participated in by TCS molecules to the total number of elementary reactions

Figure 5

Proportion of elementary reactions participated in by the SiHCl3 molecule to the total number of elementary reactions at the pyrolysis temperature between 1000 and 2000 K.

Apparently, WTCS did not monotonously change with the increase of the pyrolysis temperature, whose minimum achieved 60%, while the maximum reached 85%, elucidating the high reactivity of TCS molecules in participating in chemical reactions. It should be noticed that the presence of basically an equal number of elementary reactions involving a reactant and a generator showed that the TCS molecule could also act as one kind of extremely reactive intermediate and not just as an initial reactant during the dynamics, totally different from the traditional cognition. The distinguishing of this phenomenon should be attributed to the application of the ReaxFF MD method itself.

Analysis of the Evolution Trend of Energy Terms

In general, the overall energy in the ReaxFF system was derived from the individual contributions from a variety of partial energy terms in the distinguishing of whether it was based on BO or not. The BO-dependent terms such as bond energy, valence angle, lone pair, conjugation, and torsion energy affected the appropriate disposing of the properties of the preferred configurations of atomic and corresponding molecular orbitals, while the independent ones acted to deal with the nonbonded interactions such as van der Waals and Coulomb energy. Therefore, it was supposed that some of the hidden interactions between the distribution of intermediates and the pyrolysis temperature during the simulations were indicated by the variation of these energy terms described above. To study the validity of this inference and give a more comprehensive description of the pyrolysis kinetics of TCS simultaneously, the detailed analysis of some energy terms including bond energy, van der Waals, Coulomb energy, and torsion energy is provided in the following part.

In particular, bond energy describes the energy necessary to dissociate an existing bond and is responsible for the evaluation of bond stability. Once a certain type and number of atoms were preset in the simulated system, it could be deduced that both the class and quantity of the newly generated bonds were identical to a large extent, provided that the connectivity between atom pairs was consistent with each other. In other words, the total bond energy of the system basically remained unchanged under the condition that one and only one bond was formed in an average of two atoms. However, the total bond energy would decline supposing that more of the molecules possessing long linear structures, even polycyclic structures, were created in the reaction space resulting from the shaping of less than one bond between either of the averaged two atoms. At this point, there existed a close relationship between the total bond energy and the quantity of the bond states (bonding and breaking) in the simulated system. Just as the time evolution of the total bond energy at varied temperatures is shown in Figure a, the total bond energy in the isothermal pyrolysis system of 1000 K maintained a gradual downward tendency from the beginning to the end, which subsequently went through a gradual decline, followed by a brief ascent, and eventually a reduction in the system under 1300 K. As to the system under 2000 K, an asymptotic value of the total bond energy was reached after transient going down and going up with time evolution. Here, the solid line corresponds to the fitting curve, while the square, circle, and triangle in the open format represent the actual value, similar to those shown in Figure . These variation trends of the total bond energy illustrate that much smaller molecules preferred to aggregate into Si clusters in the system under 1000 K than that under 1300 and 2000 K, resulting from the deduction that the emergence of Si clusters in diverse sizes might make significant contributions to the continuous reduction of the total bond energy. It was not likely that Si clusters would appear in a relatively large size during the simulation under high temperatures because of the asymptotic value of the total bond energy, especially for the system above 1300 K. This means that the ascension of the isothermal pyrolysis temperature would reduce the possibility of forming Si clusters in a large size, agreeing well with the results shown in Figure and molecules shown in Figure . Clearly, the isothermal pyrolysis temperature plays a significant role in determining the distribution of the intermediates. The possibility of avoiding the yield of Si clusters in a relatively large size could be realized by the application of a low temperature in the pyrolysis system of TCS.

Figure 6

Time evolution of bond energy (a) and van der Waals energy (b) for the temperature of the SiHCl3 pyrolysis system under 1000, 1300, and 2000 K.

Figure 7

Time evolution of Coulomb energy (a) and torsion energy (b) for the temperature of the SiHCl3 pyrolysis system under 1000, 1300, and 2000 K.

In the system with a changing connectivity, ReaxFF computed the interactions including van der Waals and Coulomb interactions between every pair of atoms irrespective of the connectivity with the aim of describing the nonbonded interactions, and any excessive close-range nonbonded interaction was avoided by the inclusion of a shielding term. Broadly speaking, the van der Waals interaction accounted for the attraction that existed between molecules and was tightly associated with the average distance between molecules, the smaller value of which indicated the larger averaged distance between molecules. Similarly, the change of the van der Waals interaction with time evolution differed with the specific pyrolysis temperature in the application (Figure b). The relatively larger van der Waals interaction occurred in the system under high temperatures of 1300 and 2000 K than under a low temperature. As to the concrete variation tendency of the van der Waals interaction with the reaction progress, the slow reduction trend under 1000 K faded away with the ascend of temperature, for example, to be replaced by the one in long asymptotic trend later on a brief declare under 2000 K, the one obtained under 1300 K was somewhere between that of 1000 and 2000 K. The lower temperature below 2000 K and the existence of large amounts of TCS molecules at the beginning of the reaction under 2000 K were responsible for producing Si clusters, leading to the extension of the average distance between diverse molecules and radicals because of the application of the NVT ensemble, further directly undermining the mutual attraction of molecules to each other, which manifested as the van der Waals interaction in turn. The gradual increase to a high-temperature (2000 K) with the ongoing simulation blocked the continuous growth of these Si clusters with respect to the size and quantity, even decomposing them into small molecules or radicals with a nearly uniform size, causing the fluctuation of the total van der Waals interaction to be centered on one asymptotic value. Therefore, it follows that the essence of the consequence of the pyrolysis temperature on the van der Waals interaction is mainly implied in the formation of Si clusters. In addition, as to the system under 1000 K, the reduction of the van der Waals interaction to a slow rate with the reaction progress validated its weaker interaction than that of the bond energy in identical environments including intermediate distribution and chemical events.

In the case of another nonbonded interaction, the temperature dependence of the Coulomb energy on time evolution was nearly kept identical to that of the van der Waals interaction. However, as shown in Figure a, compared with the change of the van der Waals interaction, the main difference between the variation trends of the Coulomb energy lied in its high concentration of actual value around the fitted value. This might come from the fact that the charge and distance between molecules together had a consequence on the Coulomb energy. Also, the distance between molecules was inversely proportional to the value of the Coulomb energy. The rearrangement of the partial charge of atoms could occur due to the formation of new intermediates such as free radicals, and the rearrangement degree was significantly based on the pyrolysis temperature. Additionally, the comparison of the sharp reduction of total Coulomb energy under 1000 K with the slow reduction of total van der Waals interaction under 1000 K illustrated the stronger interaction of the Coulomb energy than the van der Waals interaction on the premise of the same reaction conditions. Thus, it was deduced that the production of the Si aggregating species in dynamic with time as the reaction progresses also made the main contribution to the variation in total Coulomb energy.

Since ReaxFF had been extended to partition the total energy of the system into individual contributions from several partial energy terms, the declaration of some energy terms might naturally cause the ascension of other energy terms to balance the total energy. Taking torsion energy as an example, as shown in Figure b, the variation trend with time evolution was entirely different from the above discussed bond energy, van der Waals, and Coulomb energy, especially for the results presented in the system under 1000 K. In a nutshell, the definition of torsion energy was rooted in characterizing the energy change associated with the distortion of the molecular skeleton resulting from the rotation of a single bond. The molecular skeleton of one asymmetrical molecule was supposed to be seriously distorted by rotating its single bonds on the condition of containing large amounts of single bonds, resulting in the jump rise of torsion energy. That is to say, the generation of several Si clusters in a relatively large size during the simulation would be responsible for the sharp increase of torsion energy since the torsion energy described in Figure b shows the aggregation of torsion energy for all of the instantaneous molecules and not the individual value for one molecule. Therefore, the non-stop occurrence of the polymerizing small molecules or radicals in Si clusters in large size with a polycyclic or linear structure significantly promoted the increase of the total torsion energy with time as the reaction progresses. For the same reason, the ascension of total torsion energy at the preliminary stage of the system under 2000 K was attributed to the generation of molecules with many atoms. It followed that the final asymptotic total torsion energy lasting for long reaction time should be ascribed to the decomposition of sophisticated molecules into small ones caused by the application of a high temperature (2000 K). This might be because the BO-dependent interactions of torsion energy were connectivity-dependent in the meantime so that the absence of energy contributions from the torsion energy took place upon bond dissociation. The presence of not-so-obvious variation tendency of total torsion energy in the system under 1300 K was derived from the gentle size change of the momentary intermediates. As a result, the change of total torsion energy with time evolution was also tightly associated with the dynamics of the distribution of the intermediates. In some sense, the product distribution in the present work is consistent with the results provided in the investigation of the gas-phase reactivity of chlorosilane,[8] where the disilane mechanism was considered as dominant at low gas temperatures while the radical mechanism was regarded as the control at high temperatures, since a substantial number of chlorosilanes were generated at low temperatures in the present simulations.

In conclusion, the variation of all the segmental energy terms discussed above was excellently in line with the distribution of instantaneous intermediates, which was greatly affected by the isothermal pyrolysis temperature in turn. To illustrate this point of view more clearly, the intermediates including molecules and radicals generated instantly during the simulation were lumped and are shown in Figure . In particular, these intermediates were classed as Si0, Si1, Si2, Si3, Si4, Si5–Si10, and Si11–Si17 clusters based on the maximum number of Si atoms in a single molecule, simultaneously considering the total quantity of these lumped species. As shown in Figure , it should be noted that the scattered open form shaped in a square, a circle, a triangle, and a star represented the real value while the solid lines indicated the resulted value fitted by Gaussian. For convenience, the entire reaction time was partitioned into three periods in order to efficiently describe the temporal interval available for generating these lumped species, namely, the first, second, and third periods finished at 300 ps, 700 ps, and the end of simulation, respectively. On the whole, the formation of Si2 clusters and the consumption of Si1 clusters in abundance dominated the first period of the simulation. It was extremely likely to produce diverse Si2 clusters from the bonding between Si1 clusters. When the reaction reached the second period, the peak yield of Si2 clusters was achieved, the Si1 clusters gradually disappeared, and Si3 and Si4 clusters were obtained in large numbers. On the basis of these phenomena, the generation of Si3 and Si4 clusters could be attributed to the reaction between the Si2 clusters and Si1 clusters, where Si2 clusters play the vital role. At the last third period, the clusters in large size including Si5–Si10 and Si11–Si17 species began to appear, along with the consumption of Si2 clusters, the relatively large generation of Si0 species, the mild variation of Si3 and Si4 clusters, and the nearly unchanged Si1 species. It seemed that the participation of Si2, Si3, and Si4 clusters was crucial and essential for the appearance of Si clusters in large size. Since quite a large number of Si atoms lying in single molecules belong to Si5–Si10 and Si11–Si17 clusters, it did not mean a whole lot in terms of the absolute number of one particular cluster. In addition, the chemical events associated with the formation of Si clusters in large size gave birth to Si0 species such as H*, Cl*, HCl, HCl2, and H2Cl2 in the meantime. However, the contributions of these Si0 species to the momentary variation of total partial energy terms were lower in comparison with the Si clusters in large size. On one hand, it followed that the distribution of intermediates along the reaction progress was actually in accordance with the change of partial energy terms. That is to say, the analysis of the variation of partial energy terms could give a brief introduction to the composition of intermediates formed at the moment. On the other hand, the desirable intermediates that put much interest in could be obtained by adjusting the total simulation time and the employed temperature.

Figure 8

Variation of the number of intermediates with time evolution under the pyrolysis temperature of 1000 K: (a) Si0, Si1, Si2, (b) Si3, Si4, Si5–Si10, Si11–Si17.

Reaction Pathway Properties

To get a better knowledge of the distribution of intermediates and not just the lumped ones presented above, the particular molecules and radicals with a relatively high frequency of occurrence that appeared in the system under 1000, 1300, 1600, and 2000 K are presented in detail in Table . Here, for convenience, the label of top was employed to denote the first n of molecules or radicals sorted in all of the intermediates according to the comparison of their total occurrence frequency during the whole simulation process. Since the intermediates presented relatively high occurrence frequency and were supposed to be equipped with representativeness in terms of the distribution of intermediates, only the representative intermediates from top1 to top10 are provided in Table . Obviously, the pyrolysis temperature played a crucial role in the total occurrence frequency of the intermediates, where the small and simple ones such as SiHCl2, H*, and HCl preferred to appear at a low temperature while the relatively sophisticated ones such as Si3H3Cl9 and Si3H3Cl8 were inclined to emerge at a high temperature. The sophistication of one molecule was not only embodied in possessing relatively more atoms but also manifested in the conflicting spatial structures under a uniform formula. Remarkably, both the high reactivity and stability of the intermediate made collective contributions to the total occurrence frequency of one particular molecule and/or radical. Therefore, the fact that the maximum Si atoms in one single intermediate from top1 to top10 only reached 3 demonstrated the unstable and short-lived properties for most of the Si clusters in large size. Presumably, the maintainability of high frequency of occurrence for SiHCl4 and Si2H2Cl6 even under poor stability could be attributed to extremely high reactivity that was continually involved in chemical events. This means that, for instance, most of the TCS molecules would first grow into the intermediate state of SiHCl4 to further react in chemical events, acting as the inevitable reaction pathway and resulting in the greatly uphill frequency of occurrence for SiHCl4 in turn. The appearance of unstable Si2H2Cl6 molecules in a high frequency of occurrence was to characterize the reversible reaction that occurred too often between two TCS molecules. As a whole, the bridge bond of Si–Cl–Si appeared in nearly all of the Si2 and Si3 clusters, sufficiently validating the essential factors bringing these lumped intermediates into existence. Although the lifetime of molecules and radicals tended to be cut short in the premise of the coverage of one single Si atom by five bonds, its reactivity went through improvement in the other way round. The occurrence of small molecules and radicals such as SiHCl2, HCl, H*, and SiCl4 was speculated to be rooted in the self-decomposition of TCS molecules ahead of its encounter with the other small molecule or radical containing a Si atom.

Table 1

Top10 Intermediates in the SiHCl3 Pyrolysis Process at Different Temperatures: 1000, 1300, 1600, and 2000 K

In order to explore the consequence of pyrolysis temperature on the distribution of intermediates on the whole, the molecules and radicals sequenced in the averaged first five among the overall temperature range (1000–2000 K) are presented in Table , namely, the species began at top1 and ended at top5. Apart from the original existing TCS molecules, Si2H2Cl6 and SiHCl4 molecules with poor stability and high reactivity were proved to be the first and foremost intermediates regardless of the employed temperature for pyrolysis. The presence of HCl molecules in these first five intermediates might be rooted in their relatively high stability and the character of not easily being captured by other molecules and/or radicals because of their small molecular volume. The cleavage of the Si–Cl bond in Si2H2Cl6 would naturally give birth to Si2H2Cl5 along with HCl molecules. Of course, this was just one possible reaction pathway to create Si2H2Cl5; there were also other available routes to make certain contributions for the complicated overall reaction network.

Table 2

Averaged Top5 Intermediates in the SiHCl3 Pyrolysis Process at Temperature between 1000 and 2000 K

To characterize the above-mentioned intermediates in quantification, the proportion of the quantity of these top5 intermediates to the total quantity of all the generated intermediates is shown in Figure and could be simplistically depicted as the following expression

Figure 9

Occurrence proportion of the top5 intermediates (scatter line in blue) to the total occurrences of all intermediates (column in dark yellow) in the SiHCl3 pyrolysis system at various temperatures.

It was postulated that these first five intermediates had great correlation with the chemical events that occurred in the pyrolysis system of TCS owing to their high proportion that fluctuated between 76 and 95%. Notably, this proportion and the total occurrence times of all intermediates gradually went downhill in correspondence to the declaration of the pyrolysis temperature. Even though the occurrence times of all the intermediates in summary was the lowest among the systems under various temperatures, this proportion reached its minimum in the system under 1000 K since there were so many kinds of molecules and radicals produced in the system other than the referred first five intermediates. The comparatively small occurrence times of all the intermediates at the system under 1000 K would naturally remind us of the continuous generation of Si clusters at the third stage of the reaction process. It is worth mentioning that the gap between this proportion and the total occurrence times of all intermediates tended to be extremely small, indicating that the consequence of further increasing the pyrolysis temperature on the primary intermediates was probably negligible until reaching one certain high temperature. On the other hand, the concentration of the chemical events involved in these mentioned first five intermediates became enhanced by rising the pyrolysis temperature to one particular value. In conclusion, these phenomena indicate that prompting the pyrolysis products centered at certain species could be improved by elevating the isothermal temperature.

As to the specific reaction pathways for one specie, the pyrolysis temperature also made contributions of significance. For convenience, only the product of SiHCl2 was focused on to proceed with the detailed analysis, even though so many intermediates were produced during the simulation. There were only three primary reaction channels relevant to SiHCl2 obtained in the system under a high temperature of 2000 K, which were provided in the following

However, the reaction disappeared in the system under 1000 K, even though substantial SiCl4 was generated in the dynamics. Besides, the reactions G1 and G2 only functioned as the primary ones but not the exclusive ones. It was, however, significant to point out that there were many chemical channels involved in SiHCl2 and chlorosilanes, such as the reaction channel depicted as followswhich could be ascribed to the fact that substantial chlorosilanes existed in the system. In light of these results, it was found that the pyrolysis temperature could determine the particular reaction pathways for the transformation of one product.

Effect of Pyrolysis Temperature on SiHCl3 and Si2H2Cl6

Based on the above discussion, molecules of TCS and Si2H2Cl6 were proved to belong to the most active intermediates, especially for TCS, prompting this part to probe into their dynamic evolution in terms of conversion and occurrence frequency. According to the variation trend depicted in Figure , not only the final number of TCS molecules at the finishing point of pyrolysis process but also the conversion of TCS molecules during the entire simulation were sensitive to the pyrolysis temperature. However, the consequence of temperature on the final number and conversion of TCS molecules was contradictory, resulting from the constant number of TCS molecules at the beginning of simulation, yet this consequence grew into virtual absence under a high temperature. Specifically, as to the TCS molecule, the final number gradually went up, whereas the conversion inversely went down with the elevation of temperature. Besides, two sudden jumps and declines relevant to the change of the final number and conversion, respectively, could just be observed obviously under the temperature lower than 1600 K, illustrating the very tiny influence imposed by the high temperature. For instance, the conversion reached its maximum of 62% under 1000 K; nevertheless, it was kept lower than 20% in the temperature range above 1500 K. Just as in the above discussion, the aggregation of small molecules and radicals into Si clusters a large size in the system under a low temperature would simultaneously destroy the TCS molecules through both the direct and not direct pathways, giving birth to the relatively small number of TCS molecules in existence at the end of simulation. Even though these Si clusters possessed poor stability, they just transformed into other Si clusters along with a negligible change of the polycyclic structure rather than decomposing into individual small molecules in one time. Therefore, the selection of a high pyrolysis temperature became essential to boost the degree of pyrolysis reaction of TCS molecules.

Figure 10

Final number and conversion of SiHCl3 at the end of the pyrolysis process with various temperatures.

In general, the frequency occurrence of Si2H2Cl6 molecules was possibly rooted in their extreme poor stability and the resulted high reactivity since five bonds encircled one unitary Si atom of Si2H2Cl6 molecules. As Figure shows, the low temperature was favorable for the generation of Si2H2Cl6 molecules, speculating from its relatively high occurrences’ proportion to total occurrences of intermediates. This might be because substantial amounts of opportunities for the encounter between two TCS molecules were available for the formation of Si2H2Cl6 molecules deriving from their primary transformation pathway of the revisable reaction between two TCS molecules. On one hand, it is postulated that TCS molecules preferred to frequently polymerize into the highly reactive Si2H2Cl6 molecules prior to their clustering under a low temperature, resulting from relatively larger amounts of free TCS molecules in the system. On the other hand, the energy for active TCS molecules to proceed in other chemical events apart from forming Si2H2Cl6 molecules might not be large enough under some given simulation condition. In practice, however, the occurrence proportion of Si2H2Cl6 molecules began to decline little by little; yet, that of TCS molecules began to upgrade by degrees along with the advance of temperature. The high temperature made a contribution to provide enough energy for promoting the chemical events associated with the pyrolysis of TCS molecules to occur step by step not merely the reversible reaction between two TCS molecules, resulting in the reduction of the occurrence proportion of Si2H2Cl6 molecules in the meantime and the presence of peculiar radicals that could not be observed under a low temperature such as SiCl2. These phenomena illustrated that even though there were many kinds of molecules and radicals during the pyrolysis process, the absolute number of each individual kind of molecule and radical was relatively small in comparison to TCS molecules. Thus, some of the reaction pathways under a high temperature were completely different from those under a low temperature. However, the accurate comparison of the imposed influence of temperature on the transformation pathways of one specific molecule or radical was invalid based on the truth that the long pathways were flexible. There was a special case that the molecule or radical of interest possibly might not be formed at a low temperature. In addition, the pyrolysis temperature was inclined to preset at the low one in order to realize the aim of acquiring more Si2H2Cl6 molecules.

Figure 11

Proportion of SiHCl3 and Si2H2Cl6 occurrences to total occurrences of all intermediates at different pyrolysis temperatures.

Concluding Remarks

To explore the distribution of intermediates as well as the consequences exerted by the pyrolysis temperature, ReaxFF was carried out to perform a series of MD simulations on SiHCl3 molecules under isothermal conditions. Our simulations demonstrated that the distribution of intermediates varied with both of the simulation period and the pyrolysis temperature do have significant influences on the essence of pyrolysis products. In general, complex intermediates containing more Si atoms preferred to appear under a temperature lower than 1300 K such as Si clusters with a polycyclic structure, and the simple ones tended to emerge under a high temperature, such as SiHCl2, H*, HCl, and so on. For instance, as to the maximum number of Si atoms in a single molecule, two jump rises were observed under a temperature lower than 1300 K. Particularly, the appearance of lumped Si clusters at the later period of the simulation was drawn by a low temperature, which served likely as nucleation points for further aggregation of small molecules and/or radicals until the formation of a polycyclic structure with a large size and a short life. The characters of activity and lifetime work together to determine the occurrence frequency of one individual molecule and radical. These findings were in good agreement with the observations of partial energy terms such as the BO-dependent ones of bond energy and torsion energy as well as the BO independent ones of van der Waals and Coulomb energy, indicating that the variation of segmental energy terms could be attributed to the dynamic evolution of intermediates in terms of the definition of these energy terms. In addition, the conversion of SiHCl3 molecules was also sensitive for pyrolysis temperature, while their consumption with time evolution presented a linear fitting only under a temperature lower than 1300 K. Furthermore, the first 10 ones sequenced in all kinds of the intermediates that formed during the whole simulation process were described in detail in the present work, especially under the temperatures of 1000, 1300, 1600, and 2000 K. The results presented here elucidated that the reaction environment was more complicated under a low temperature than that under a high temperature. Interestingly, the SiHCl3 molecule could come into action as an active product not just the initial reactant. Our future works will focus on illustrating the micro-reaction mechanism of the integrated CVD process containing SiHCl3, H2, and silicon substrate.

In general, the above simulation results in discussion revealed that the H2 environment and/or silicon substrate was the indispensable reaction conditions for the generation of monatomic silicon. Besides, the acquirement of the influence of reaction temperature on the dynamic evolution of the main intermediate distribution provided some valuable suggestions concerning the reaction temperature applied in the ReaxFF MD simulation of the CVD process of Si from SiHCl3. Also, the announcement of the intrinsic transformation of the Si element in SiHCl3 before its deposition on the silicon substrate provided help for distinguishing the role of a H2 environment in the gas-phase reaction of the CVD process in precise. Therefore, this work provides a guidance for investigating the micro-reaction mechanisms of the CVD process of high-purity polysilicon.

Calculation Methods

ReaxFF Molecular Dynamics

In general, the ReaxFF practically preserves the accuracy of QM and allows MD for computational costs as low as for simple force field at the same time, building one bridge that connected QM and MD methods, enabling it to calculate large and complex reaction systems with thousands of atoms or more for periods of nanoseconds or longer. Differing from the traditional nonreactive force field, the BO-dependent ReaxFF can supply precise descriptions of bond dissociation and creation continuously without setting the reactive sites and reaction pathways in advance. As a function of interatomic distances, BOs are calculated in every iteration step[17] since it determines the connectivity between atom pairs, playing a role as the fundamental difference between the traditional nonreactive force field and ReaxFF, allowing for bonds to cleave and form during the dynamics. In other words, nonfixed connectivity assignments for the chemical bonds are applied in the ReaxFF. A detailed description of the Si/C/Cl/H ReaxFF parameters employed in the present work is provided in the Supporting Information, deriving from the successful calculation in the adiabatic reactive dynamics of silicon carbide growth from a methyltrichlorosilane precursor.[18]

ReaxFF divides the overall system energy into the contributions expressed as follows[19]

Clearly, the overall system energy is the sum of all the partial energy terms, including bond Ebond, three-body angle Eval, four-body torsional angle Etor energies, and non-bonded energies, such as van der Waals Evdw, Coulomb Ecoul interactions, as well as the hydrogen bond EHB contribution. In addition, the terms of Elp, Eover, Eunder, Epen, Ecoa, and Econj are the correcting energy, namely, Elp for the appearance of lone pair, Eover and Eunder for over and under coordination of atoms associated with their valency, penalty energy Epen for stabilizing one three-body system in the case of the creation of two double bonds between the centered atoms, Ecoa and Econj, for the presence of conjugated chemical bonds.

BO Cutoff

In order to understand the chemical events in terms of molecules relevant to the pyrolysis process of TCS, the critical point lies in how to recognize one as a molecule, which in reverse is dependent on the identification of the bond, namely, the connectivity between any two atoms. In general, two steps are essential to compute BOs,[16] which determine the connectivity of each atom pair, including the first calculation of the original BO (BO′) and the subsequent modification of BO′ by the parameters of quantum chemistry regression. Unlike BO, BO′ can be directly computed from instantaneous interatomic distance r as described in the following eq where BOij′σ, BOij′π, and BOij′ππ represent the BO of single, double, and triple bonds, respectively. Additionally, r0σ, r0π, and r0ππ mean the balance distance of single, double, and triple bonds, respectively. When two atoms are closer than the specific cutoff distance, they are supposed to belong to the same molecule. Based on this, the absolute threshold of BO is necessary to judge the binding state between two atoms, which can be presented as the BO cutoff varying with the elements involved. As to one atom pair’s BO, a greater value than the BO cutoff indicates the formation of a bond between this atom pair, and the smaller one is on behalf of the nonbonded state. In the present work, the distribution of BO during the simulation is plotted in Figure , resulting from the update of BO in each time step.

Figure 12

Bond order distribution of chemical bonds: (a) Si–Si, (b) Si–Cl, (c) Si–H, (d) H–Cl, (e) H–H, (f) Cl–Cl. The red solid and blue dashed lines represent fitting curve and the corresponding bond order cutoff, respectively.

Based on the expression of free energy ΔG = –RT ln P(x), the P(BO) minimum of the BO distribution curve (Figure ) corresponds to the maximum of ΔG and can be employed as the BO cutoffs of chemical bonds,[20] where P(BO) represents the probability distribution of BO during the overall evolution of the simulation system. Generally speaking, P(x) in the description of ΔG = −RT ln P(x) derives from the equivalent transformation of the partition function (Z) employed in statistical mechanics and the canonical ensemble, where Gibbs free energy can be described as G = −kBT ln Z, with kB and T representing Boltzmann constant and temperature, respectively. Based on the simplified physical meaning of regarding the partition function as the normalization constant of probability in numerical terms, the probability distribution of BO in the simulation is basically in agreement with the partition function, resulting in the application of P(BO) in the expression of ΔG = −RT ln P(x). Clearly, compared with other chemical bonds, the BO distributions of Si–Cl and Si–H chemical bonds are smoother, caused by more Si–Cl and Si–H chemical bonds in number in the pyrolysis system. Based on the BO variations, fitting curves in a red solid line could be plotted to give the specific BO cutoffs in blue dashed line, summarized in Table . Analysis of intermediates generated during the simulation was carried out with the BO cutoffs in Table for the identification of intermediates. Even though the selection of BO cutoffs did not play a role in the simulation itself but only the determination in accordance with chemical components, the capture of unsuccessful reactions generating very short-lived species could be ignored by the reasonably set BO cutoff.

Table 3

Bond Order Cutoffs of Chemical Bonds

chemical bond	Si–Si	Si–Cl	Si–H	H–Cl	H–H	Cl–Cl
BO cutoff	0.30	0.55	0.60	0.38	0.55	0.30

Simulation Details

To carry out TCS pyrolysis simulations, a periodic system including 100 TCS molecules was established. The configuration of the system containing the random distribution of TCS molecules at a density of 0.001 g/cm3 in a cubic box with an edge length of 282.3 Å is shown in Figure . All of the simulations were performed with a constant number of N atoms in a constant volume V where the control of temperature T was realized by using a Berendsen thermostat[21] with a damping constant of 0.02 ps. That is to say, the conditions of the canonical ensemble (NVT) were applied during the dynamics with a fixed timestep of 0.2 fs. The simulation procedure began with the energy minimization of each reaction system under the employed ReaxFF to acquire a reasonable configuration of the initial reaction model. Next, the system was simulated at a low temperature of 300 K with 50 ps for relaxing the atomic positions, followed by the heat-up simulations from 300 K to the desired pyrolysis temperature at a heating rate of 10 K/ps. Subsequently, the system was simulated under isothermal reactions with the duration as long as 1000 ps. To evaluate the temperature dependence of the pyrolysis kinetics of TCS molecules, a series set of simulations were performed at varied pyrolysis temperatures (from 1000 K in compression to 2000 K in expansion). It should be declared that the instantaneous intermediates would persistently remain in the system to participate in the further chemical events during the entire simulation. To acquire the detailed chemical reaction information as much as possible, all of the relevant data including the movement trajectory of atoms was recorded in each 0.2 ps. In addition, each system was carried out in three parallel simulations with a unique initial configuration to obtain a statistic distribution of the observed reactions and intermediates.

Figure 13

Schematic diagram of the reaction model for the SiHCl3 pyrolysis system.

All of the simulations in the present work were performed with the LAMMPS code. For greater convenience and effectiveness, VARxMD (Visualization and Analysis of Reactive Molecular Dynamics) software[22] was employed as the molecule recognition method to analyze the simulations because of its previous successful usage in reading the chemical events in the extremely complex pyrolysis system of coal[13,23,24] and biomass[25,26] simulated by the ReaxFF MD method.