Article pubs.acs.org/IECR
Extending Reaction Mechanism Generator to Silicon Hydride Chemistry Belinda L. Slakman,† Harsono Simka,‡ Harinath Reddy,¶ and Richard H. West*,† †
Department of Chemical Engineering, Northeastern University, Boston, Massachusetts 02115, United States TCAD, Logic Technology Development, Intel Corporation, Santa Clara, California 95054, United States ¶ TCAD, Logic Technology Development, Intel Corporation, Hillsboro, Oregon 97124, United States ‡
S Supporting Information *
ABSTRACT: Understanding the gas-phase chemistry of silicon hydrides is the first step to building a realistic kinetic model for chemical vapor deposition (CVD). Functionality for thermodynamic and kinetic data estimation of silicon hydrides was added to the open-source software Reaction Mechanism Generator (RMG). Using the updated RMG, a detailed kinetic model was built for SiH4 thermal decomposition. The generated model was used to perform reactor simulations at various process conditions for comparison to prior SiH4 decomposition experiments in a flow tube. Results show that the RMG-generated model can reasonably replicate experimental results for SiH4 concentration profiles at different temperatures and residence times. While the effect of changing initial SiH4 concentration is not captured, a first pass sensitivity analysis reveals that reasonable errors in reaction rates could contribute to the discrepancy.
■
INTRODUCTION The gas-phase chemistry of the precursors to chemical vapor deposition (CVD) directly affects the yield and quality of the solid materials produced. Silane (SiH4) CVD has been wellstudied experimentally and theoretically,1−5 but all reaction pathways, including those involving radical species, should be considered in modeling and prediction of the CVD process. Complete, detailed models of the reactions of these precursor gases can provide necessary insights into the semiconductor industry, and studying SiH4 can serve as a proof-of-concept that can be extended to other novel gas precursors. Detailed kinetic models also allow predictions to be made at different operating conditions. Furthermore, a framework for studying gas-phase silicon hydrides provides the first step in studying surface reactions involved in CVD. Mechanisms for SiH4 thermal decomposition have been generated from experimental data and by hand,6−8 with the largest manual model including about 220 chemical species and 2600 chemical reactions.9 However, building these large mechanisms by hand is tedious, and automatic mechanism generators can be used to build the models faster and minimize errors. Wong et al.10 previously used automatic mechanism generation to build a model for silicon nanoparticle formation from SiH4 decomposition. In that study, automatic mechanism generation software developed by the Broadbelt group, which © 2016 American Chemical Society
uses a rate-based termination criteria to limit model size, was employed. Reaction types included silylene (SiH2) insertion into Si−H and H−H bonds (and its reverse, SiH2 elimination); silylene to silene isomerization; and ring opening and closing isomerizations. The study included a comprehensive analysis of the different factors contributing to silicon particle clustering. However, radical pathways were not included in the built models. Several decades ago, Purnell and Walsh suggested that decomposition of SiH4 to SiH3 and H plays a role despite not being the most dominant pathway.6 An earlier study by Emeléus and Reid11 suggested the role of SiH3 becomes more important when Si2H6 and Si3H8 are used as precursors. The involvement of radicals may also be more important at certain CVD conditions,12,13 for example, in low temperature CVD in silane plasma.14 SiH3 was found to be the dominant radical species at steady state in silane plasmas in a study by Robertson et al. 15 Watanabe et al. later demonstrated that the concentration profiles of Si, SiH, and SiH2 radicals greatly affect particle growth in these systems.16 Received: Revised: Accepted: Published: 12507
June 22, 2016 October 31, 2016 November 1, 2016 November 1, 2016 DOI: 10.1021/acs.iecr.6b02402 Ind. Eng. Chem. Res. 2016, 55, 12507−12515
Article
Industrial & Engineering Chemistry Research
database of thermodynamic and kinetic parameters. Radical reaction types, which already existed in RMG but lacked data for silicon, and new reaction types specifically important for silicon hydride chemistry, can both be proposed with the updated RMG. Using the new data in RMG, a model for SiH4 thermal decomposition was built. The resulting model was used to simulate a flow reactor, and these simulations were compared to SiH4 flow tube experimental data obtained from Onischuk et al.,22 including data for SiH4 and Si2H6 concentration profiles with time, and at different temperatures, residence times, and initial concentrations of SiH4.
The Reaction Mechanism Generator (RMG) is an opensource, automatic mechanism generation software developed at MIT and Northeastern University.17,18 RMG has typically been used to build models involving gas-phase hydrocarbon combustion and has functionality for chemical species containing carbon, hydrogen, oxygen, sulfur, and nitrogen. RMG stores chemical species as graphs, with atoms as nodes and bonds as edges, and identifies matching species using graph theory methods.19 To generate reactions, RMG uses reaction families, each representing an elementary chemical reaction type. Each family consists of a template, which defines the reacting molecular group(s) that can participate in that family, and a recipe, which defines a set of actions occurring in the reaction (e.g., bonds forming or breaking, radical electron count increasing or decreasing). An example of the template and recipe for the hydrogen abstraction reaction family is shown in Figure 1. To
■
UPDATING RMG’S DATABASE To enable RMG to generate detailed kinetic models for silicon hydrides, small updates to the RMG source code were made, and RMG’s database was updated with both kinetics and thermodynamics data. New kinetics can be added to RMG’s database in a reaction library and/or as training reactions in a reaction family. Rates in a reaction library are used, without modification, whenever the reactant (or product) species are found in the model core. This is useful for reactions that do not behave nicely as part of a reaction family, e.g., small molecules that go against the trend or reactions that do not match an existing recipe, and for pressure-dependent reaction rate expressions that should bypass RMG’s Master Equation code. On the other hand, reaction kinetics added as training reactions to a reaction family are placed into the hierarchical tree of estimate rules, using the most specific functional group definitions possible for those reactants. They are then used when completing the more general nodes in the tree by averaging. Training reaction rates must provide high-pressure limit kinetics for an elementary reaction of a specific reaction family, and will influence the estimates of other similar reactions generated by the reaction family. This is the preferred way to add kinetic data to RMG’s rule-based estimation database.17 Published reaction rates have been added to reaction libraries, which are preferentially used during RMG simulations. Additional rates were calculated using transition state theory. Some published and calculated rates have been added to the training data for four reaction families and will be used to estimate reaction rates in these families where the exact rates are unavailable. Published, calculated, and group additive thermodynamic data for silicon hydrides have also been added to RMG’s database. Kinetics Data. In this work, two reaction rate libraries were added to RMG’s database. The libraries are based on one experimental study of SiH4 CVD2 and one theoretical study involving pressure-dependent rate calculations for mainly Si2 species.5 These reactions and their rate parameters can be found in the Supporting Information. Two reaction families were added to RMG’s database: silylene insertion, in which a silicon atom with a lone pair can insert into a silicon−hydrogen or hydrogen−hydrogen bond, and silylene-to-silene isomerization, in which a hydrogen atom migrates from an sp3 bonded silicon atom to a silicon atom with a lone pair, and a double bond is formed. These families are based on the reaction types used in a previous work on model generation.10 The kinetics data in these reaction families come from rates calculated by Adamczyk et al.23−25 using G3// B3LYP.26
Figure 1. Example of the template and recipe for the hydrogen abstraction reaction family. In the recipe, starred atoms participate in each action; “S” refers to a single bond being formed or broken, and “1” indicates that the radical count increases or decreases by 1.
determine the rate of a reaction within a family, its reacting functional groups are identified. The functional groups are organized into a hierarchical, molecular structure group tree for each family. The trees are descended to find the most specific description of the reaction’s reacting functional groups, and a database of rate rules is searched. If an exact match for the specific groups is found, those kinetics will be used; if not, RMG’s algorithm will use the data from more general molecular structure groups. If there are still no data found, the data from the “children” of the matching generalized groups will be averaged and used as an estimate for the reaction rate parameters. A detailed description of the algorithm is provided elsewhere.17 For thermodynamic parameters, Benson’s group additivity method20 is used when there are no experimental data or data from high level calculations available in the database. RMG uses a “core/edge” reaction model to determine which reactions are important enough to include in the model.21 When new species are generated, they are put on the edge. If, during a simulation at conditions specified by the user, the rate of formation of an edge species is high enough according to a user-defined tolerance and a characteristic rate, that species is brought into the core, and all the reactions between it and the other core species are generated. The characteristic rate depends on the rates of formation of the species in the current core. The core is then iteratively expanded until a termination criteria is reached, based upon either a target time or an initial species conversion. In this work, a new element, silicon, was added to RMG. Specifically, data for silicon hydrides were added to RMG’s 12508
DOI: 10.1021/acs.iecr.6b02402 Ind. Eng. Chem. Res. 2016, 55, 12507−12515
Article
Industrial & Engineering Chemistry Research Table 1. Reaction Families Used to Generate Mechanisms for Silicon Hydrides in RMG
The third reaction family considered in previous mechanism generation for silicon hydrides was ring closing reactions (and its reverse, ring opening), in which a silylene molecule with at least three silicon atoms isomerizes to form a ring structure. Two of these reactions (for Si3H6 and Si4H8) were added into a reaction library. Because the current work does not consider larger molecules, we did not create an RMG reaction family for arbitrarily large rings. In addition, the training reaction databases of two existing reaction families in RMG, hydrogen abstraction and radical recombination, were updated with silicon hydride-containing reactions and kinetic parameters,13,27 found in the NIST kinetics database.28 For hydrogen abstraction reactions, ten additional rates were calculated. The reactant and transition state geometries were optimized using M06-2X29,30/6-311+G(3d2f),31−35 with a clear saddle point found for each reaction. The CanTherm package, which applies conventional transition state theory, was used to determine the Arrhenius kinetic parameters.36 Reliability of these rates were tested by comparing to available published experimental and theoretical rate calculations for hydrogen abstraction reactions of silicon hydrides. The templates for these four reaction families, which are used to generate the SiH4 thermal decomposition models, are given in Table 1. Since these general reactions are reversible, reactions based on the reverse templates (i.e., silylene elimination) are also included in the model. Thermodynamics Data. Group additivity values for the thermodynamics of silicon hydride species (silanes, silenes, and silylenes) were previously calculated by Swihart and Girshick.9 The group values were later improved upon by Wong et al. by fitting to G3//B3LYP calculations.37 In this work, the Wong et al. values were added to RMG’s database for use in the existing group additivity scheme for calculating thermodynamic parameters in RMG. For radical species, group additivity values have not previously been determined. We used G3//B3LYP quantum chemistry calculations to generate hydrogen bond increment (HBI) values for silicon hydride radicals. HBI is a method used to calculate the thermodynamic parameters of a radical species R* from its parent molecule RH, which is the chemical species created by saturating the radical with hydrogen atoms.38 The thermochemistry is calculated according to the following three equations:
where R* is the radical species, RH is the parent molecule, and HBI are values for each of enthalpy, heat capacity, and entropy that are necessary for determining thermochemistry. Both the G3//B3LYP calculations and group additivity values using HBI were benchmarked against higher level calculations by Katzer et al.39
■
RMG MODEL GENERATION The conditions for the RMG simulation were chosen to closely mirror those of the experimental conditions in the flow tube experiment,22 with the base case conditions of T = 913 K, P = 39 kPa, and y0(SiH4) = 0.00016 in an argon bath gas. Temperature, pressure, and initial SiH4 mole fraction were varied around these values to generate a comprehensive mechanism that could be used for reactor simulations at a variety of conditions. Pressure dependence was included in some of the built models to test its effects on the model accuracy. The pressure dependence scheme in RMG has been described previously.40
■
REACTOR MODELING Once the models were built in RMG, they were tested for validity by comparison to experimental results using a constant pressure reaction model in Cantera.41 We simulated the same conditions for temperature, pressure, initial SiH4 mole fraction, and residence time given in the experiment by Onischuk et al.22 The rate constants and thermodynamic parameters used in the simulation come from the RMG-generated mechanism. Different simulations as well as sensitivity analyses were performed to understand the effects of temperature, pressure, initial SiH4 mole fraction, residence time, and model size on the SiH4 concentration profile.
■
RESULTS The thermodynamic and kinetic data calculated were compared to published data in order to validate their use in RMG’s database. After updating the database, the RMG generated models were used in simulations and compared with experimental results. The effects of changing various process conditions were considered. Kinetics of Hydrogen Abstraction Reactions. We calculated ten hydrogen abstraction reaction rates for use in RMG’s database, but only four prior published rates were available for comparison. The rate coefficients at 300 and 1000 K and Arrhenius parameters are displayed in Table 2, along with published data where available (reactions 1, 2, 3, and 5). In all cases, the rate coefficients at 1000 K are within 1−2 orders of magnitude of the published data. However, at 300 K, the discrepancy is much greater for reactions 3 and 5. Because we are considering temperatures between 800 and 1000 K for this study, the differences at low temperatures are not
° (R*) = HBI(Δf H298 ° ) + Δf H298 ° (RH) Δf H298 ° (H ) − Δf H298 Cp°(R*) = HBI(Cp°) + Cp°(RH) ° (R*) = HBI(S298 ° ) + S298 ° (RH) S298 12509
DOI: 10.1021/acs.iecr.6b02402 Ind. Eng. Chem. Res. 2016, 55, 12507−12515
Article
Industrial & Engineering Chemistry Research Table 2. Hydrogen Abstraction Rates Calculated from M062X/6-311+(3d2f) and Transition State Theory Using Cantherma no. 1
reaction
log A
Ea
log k300K
SiH4 + H → SiH3 + H2
14.8
10.7
13.1
13.9
11.7
11.9 11.2
log k1000K
source
14.2
this work
9.51
ref 42 ref 27 ref 43 (QI) ref 43 (VTST) ref 44 ref 45 this work
13.6 13.2 13.1 13.2
2
SiH3 + H2 → SiH4 + H
3
Si2H6 + H → Si2H5 + H2
4 5
6 7 8 9 10
Si2H5 + H2 → Si2H6 + H Si2H6 + SiH3 → Si2H5 + SiH4 Si2H5 + SiH4 → Si2H6 + SiH3 Si2H5 + H → SiH3SiH + H2 SiH3SiH + H2 → Si2H5 + H Si2H4 + H → Si2H3 + H2 Si2H3 + H2 → Si2H4 + H
13.4
72.3
11.3 1.03
12.4 15.9
61.6 3.55
1.65 15.5
15.7
ref 42 this work
13.8
11.1
14.0 10.1
ref 46 ref 27 this work
14.2
77.0
11.9 11.9 1.01
15.1
12.7
13.2
14.3
this work
14.8
24.6
9.6 10.9
13.4
ref 12 this work
14.5
4.55
13.9
14.2
this work
14.1
88.1
−0.97
9.45
this work
14.7
7.95
13.5
14.3
this work
13.6
126
−8.11
6.93
this work
Figure 2. G3//B3LYP calculations of ΔfH298 ° (this work) compared to high level calculations by Katzer et al.,39 for silicon hydride species with up to three silicon atoms. Structures shown for species with more than 5 kcal/mol discrepancy.
the thermochemistry for most species, which justifies our use of it for other calculations in this work. Hydrogen bond increment (HBI) values derived from the G3//B3LYP calculations of radical species and their parent molecules are given in Table 3. ΔfH298 ° for 16 silicon hydride Table 3. Hydrogen Bond Increment (HBI) Corrections Calculated with G3//B3LYPa
a
Rates were compared with those obtained from literature, where available. Units in kJ, mol, cm3, and s. Logarithms are base 10.
particularly concerning. For reactions such as reaction 1, where many experimental and theoretical data are available and are in agreement with one another, it is practical to use this agreedupon rate in RMG’s database instead of the less accurate rate calculated here. However, the comparison of these four calculated rates to the literature data shows that when data are unavailable or scarce, such as reaction 4, these DFT estimates are reasonable enough for the purposes of automatic mechanism generation, particularly at CVD relevant temperatures. If the RMG generated model were to be simulated at process conditions at which radical chemistry becomes more important, more accurate calculations should be done. Calculated Thermodynamic Data. The parity plot shown in Figure 2 reveals that calculations of ΔfH°298 using G3// B3LYP for both closed-shell and radical species compare well with the high level calculations of Katzer et al.39 Most discrepancies are less than 5 kcal/mol. Larger differences between the two methods are seen for some Si3 and multiradical species. Katzer et al. report that for proper treatment of these species, which have both multiple radicals and divalent silicon atoms, a multiconfiguration reference wave function with a postself-consistent field calculation method must be used.39 Therefore, we know that our treatment with G3//B3LYP will not be exact. However, G3//B3LYP provides a reasonable and computationally less expensive estimation of
a These corrections account for the effect of losing 1 or 2 hydrogen atoms on the enthalpy and entropy.
radicals were then calculated via group additivity and compared to the Katzer et al. values,39 illustrated in Figure 3. In this particular comparison, the group additivity values of Wong et al.37 were used to calculate the thermodynamics for the parent molecules, for the sake of consistency across species, but during an actual RMG simulation, exact species thermodynamics would be used for the parent molecules if available in the database. Because of this, Figure 3 is a reflection of both the accuracy of Wong’s group additivity values and our calculated HBI values. While most values compare reasonably well, there are again significant differences for three multiradical species which contain either double bonds or divalent silicons, for the reasons discussed above. It may be preferable, therefore, to calculate the thermodynamics of these multiradical species at a higher level and to put them in a thermodynamic library in 12510
DOI: 10.1021/acs.iecr.6b02402 Ind. Eng. Chem. Res. 2016, 55, 12507−12515
Article
Industrial & Engineering Chemistry Research
Effect of Temperature. The full pressure dependent model was further simulated at different temperatures. Figure 4b displays the plot from experiment with the pressure dependent model simulated at 863, 893, 913 (the experimental temperature), and 963 K. The simulations at 863 and 963 K show that temperatures within a ±50 K range can affect the SiH4 concentration profile enormously. Simulating the reactor at 893 K provides a close comparison to the experimental data reported at 913 K, illustrating that an uncertainty of 20 K can fully explain the difference between the model and the experiments. The main initial decomposition reaction, SiH4 ⇌ SiH2 + H2, has an activation energy of about 55 kcal/mol, with many reported theoretical and experimental determinations covering a range of about ±5 kcal/mol. A 20 K change from 913 to 893 K corresponds to only a 1.2 kcal/mol decrease in activation energy for a global reaction rate following Arrhenius kinetics with a 55 kcal/mol activation energy. In other words, the discrepancy lies well within the expected uncertainty in the activation energies. Residence Time Variation with Temperature. In the SiH4 experiment, the reactor residence time was varied by changing the volumetric flow rate of the feed gas.22 Our simulations similarly varied the volumetric flow rate into the reactor. Figure 5 compares the SiH4 concentration profile versus temperature
Figure 3. Group additivity calculations of ΔfH°298 from RMG-Py, derived from Wong et al. GAV37 and HBI corrections (this work), compared to high level calculations by Katzer et al.,39 for 16 silicon hydride radical species. Structures shown for species with more than 5 kcal/mol discrepancy.
RMG, which are used preferentially over group additivity estimates. RMG Generated Mechanisms. Two RMG mechanisms were used for most of the reactor simulations, the difference being inclusion of pressure dependent rate parameters. Tolerance was adjusted to make the size of the models roughly equal and to ensure that radical chemistry was included in both models. The mechanism including pressure dependence contains 63 silicon hydride species and 1298 reactions and the mechanism without pressure dependence contains 57 species and 578 reactions. A third (pressure-dependent) mechanism was also generated to incorporate important species and reactions at a variety of initial SiH4 concentrations for comparison to experiment. This mechanism was larger, with 83 species and 2708 reactions. Effect of Pressure Dependence. A constant-pressure reactor model was simulated with both the pressure dependent and pressure independent RMG mechanisms. The result, shown in comparison to the prior experimental results, is displayed in Figure 4a. Since the RMG model containing pressure dependence is shown to more accurately replicate the experimental result, pressure dependence was used for the remaining analysis. This result is expected, as including pressure dependent networks in the model is important at the low pressures typically used for SiH4 CVD.
Figure 5. SiH4 concentration vs temperature at different residence times, from Onischuk et al.22 (points) and from our reactor simulations.
Figure 4. Simulation results (this work) compared with SiH4 thermal decomposition experiment (data points).22 (a) Comparison of RMG models with and without pressure dependence to experimental results. Conditions are T = 913 K, P = 39 kPa, and y0,SiH4 = 1.6 × 104 (b) Comparison of RMG model, including pressure dependence, at different temperatures, compared to the experimental results at 913 K. 12511
DOI: 10.1021/acs.iecr.6b02402 Ind. Eng. Chem. Res. 2016, 55, 12507−12515
Article
Industrial & Engineering Chemistry Research
Figure 6. Simulation results for full pressure dependent mechanisms generated by RMG, compared with a mechanism generated without radical reaction families allowed. P = 39 kPa, y0,SiH4 = 1.6 × 104, and T = 913 (a) and 613 K (b).
at four different residence times. In the compared experiment, reported inlet volumetric flow rates and residence times at 843 K were not consistent given a constant volume reactor, leading to slight uncertainties in the residence time used in the simulation. The error bars in the figure represent this discrepancy only, with more information included as Supporting Information. The reactor simulations compare qualitatively well with the experimental data. In the previous section, it was noted that a 20 K difference in temperature, which is well within the uncertainty of the reaction rates, would better capture the SiH4 profile. One can imagine that the simulation data in Figure 5 would, therefore, become closer to the experimental results if the curves were shifted by +20 K; the largest discrepancy is about 40 K. Effect of Radical Reaction Families. To investigate the effect of excluding radical reaction families, a new mechanism was generated with the radical reaction families (hydrogen abstraction, radical recombination) disabled during the mechanism generation. The results for the mechanisms simulated at 913 K are displayed in Figure 6a, along with the experimental data. At these conditions, the removal of radical reactions has no noticeable effect on overall SiH4 decomposition rate. Because radical pathways are thought to be relatively more important at lower temperatures of CVD conditions, RMG mechanisms were generated at 613 K, both with all reactions included and with radical reaction families removed, and used in simulations. At this lower temperature and far slower decomposition, as shown in Figure 6b, there is still hardly a difference between the full mechanism and the mechanisms with radical reactions removed. Comparing the Si fluxes at 6 × 104 seconds, which was chosen due to the slight acceleration of the nonradical pathway at that time, we do not see a difference in the significant pathways to SiH 4 decomposition and only negligible participation by radicals. The flux diagram for the full mechanism is shown in Figure 7. SiH4 and Si2H6 concentration profiles. Figure 8 displays the concentration of SiH4 and Si2H6 with time at different initial SiH4 concentrations at 873 K. These profiles appear not to vary with the initial SiH4 concentration (8b), although the experimental results show a clear difference (8a). A sensitivity analysis was performed to assess which reactions qualitatively affected the concentration profiles of SiH4 and Si2H6. Three reactions were identified which, when the rate constants were modified by no more than 2 orders of
Figure 7. Flux diagram for Si at 6 × 104 seconds for full, pressure dependent mechanism generated by RMG and simulated at 613 K.
magnitude (102), produced noticeable changes in these concentration profiles. These reactions are as follows: SiH 2 + SiH 2 ⇌ Si 2H4
H 2 + SiH3SiH ⇌ Si 2H6 H 2 + SiH3SiH ⇌ SiH 2 + SiH4
Changing these few rates at the same time can produce concentration profiles that are close to the experimental results, as shown in Figure 8c. While the authors do not recommend changing the rates on the basis of the current evidence, the result is shown to establish that this initial built mechanism is plausible given the uncertainties. For example, the first altered reaction rate, for SiH2 + SiH2 ↔ Si2H4, was estimated by Giunta et al. to have the same rate as SiH2 + Si2H6 ↔ Si3H8, which was measured by Jasinski and Chu2,47 (this is the rate used in the RMG mechanism). However, Dollet and de Persis reported their calculated rate for this reaction to be an order of 12512
DOI: 10.1021/acs.iecr.6b02402 Ind. Eng. Chem. Res. 2016, 55, 12507−12515
Article
Industrial & Engineering Chemistry Research
Figure 8. Variation in concentration profiles of SiH4 and Si2H6 with initial SiH4 concentration at 873 K.
magnitude different.5 Therefore, it reasonable to assume considerable uncertainty in this rate. For a final model, however, modifying rates should be done carefully, with additional evidence and more knowledge about the fidelity of all rates.
conclude radical reaction pathways are much slower compared to the main decomposition pathway to SiH2 and H2 and the subsequent insertion and elimination reactions that occur. However, the new ability to automatically build more complex mechanisms for silicon hydrides in general can provide understanding which helps in selection of process conditions and/or different precursors for CVD which reduce the amount of radicals formed, thus reducing the risk for unwanted particle formation and growth in the CVD reactor. Generation of two detailed kinetic models in RMG, one with pressure dependent reactions and the other without, showed that the pressure dependent model more accurately compared to experimental data. Furthermore, using the pressure dependent model, temperature had a large effect on the SiH4 concentration profile with time. A 20 K change in temperature, which is comparable to a change of only about 1.2 kcal/mol in activation energy, brought the simulation results closer to the experimental SiH4 concentration profile. Simulations using the RMG generated model were able to reasonably replicate the effect of residence time on the SiH4 concentration profile at different temperatures. If uncertainty in temperature by about 20 K is taken into account, the simulations using our model more closely match the experiment. Comparing the SiH4 and Si2H6 profiles obtained by varying initial SiH4 concentration, the model does not match well with
■
DISCUSSION Thermodynamic and kinetic data for silicon hydrides were generated to update RMG’s database, and the computational methods used were benchmarked with available experimental and high level theoretical data. From both of these studies, it was observed that, where high accuracy data are available, they should be put into a thermodynamic or kinetic library in RMG for greater model accuracy. However, if data are unavailable, our calculations provide reasonable enough estimates for initial construction of a model. These calculations were used in RMG’s database for the hydrogen abstraction reaction family and to generate hydrogen bond increment values for calculating thermodynamics of radical species. The additions to RMG’s database allowed us to build kinetic models for SiH4 thermal decomposition. The inclusion of radical reaction families was found to be unimportant to overall SiH4 decomposition rate at both 613 and 913 K. Thus, our inclusion of these pathways does not mechanistically change the understanding of thermal decomposition of SiH4, instead corroborating the prior studies that 12513
DOI: 10.1021/acs.iecr.6b02402 Ind. Eng. Chem. Res. 2016, 55, 12507−12515
Article
Industrial & Engineering Chemistry Research the experimental data. Specifically, our simulations show no change in these profiles at different initial concentrations, whereas the experiment shows a clear dependence. We were able to obtain the qualitative trends seen in the experimental results by changing three reaction rates by amounts likely within the uncertainty of these rates; however, this approach to modifying rates is not advised, as the causes of the discrepancy may lie elsewhere. Instead, we recommend more careful calcuation of those reaction rates which can feasibly be calculated. A global sensitivity analysis could provide additional information about which rates significantly affect the SiH4 and Si2H6 concentration profiles. Further investigation by detailed measurements and calculations of these important rates could improve the accuracy of RMG’s database for silicon hydrides.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected] Phone: +1.617.373.5163. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors thank Karson Knutson from Intel Corporation for helpful discussion; RMG developers at Northeastern University and the Massachusetts Institute of Technology for help debugging; and Intel Corporation and the Department of Chemical Engineering at Northeastern University for funding.
■
■
CONCLUSION We have demonstrated a framework for extending the RMG software package to a new class of elements and chemical reactions, which allows generation of detailed kinetic models for silicon hydride decomposition. This work builds on previous efforts to automatically generate these mechanisms, making use of published thermodynamic and kinetic data. Additional data were calculated to allow radical reaction pathways to be enabled, including calculations of rates for hydrogen abstraction reactions of silicon hydrides, and hydrogen bond increment values for radical species. For silane thermal decomposition, simulations using the RMG-generated model reasonably compare to experimental results. Many other reaction families exist in RMG which can be extended to silicon hydrides, but availability of experimental or high level theoretical calculations of their reaction rates remains a challenge. This work represents an important first step in extending RMG to model chemical vapor deposition by more accurately representing the detailed gas phase chemistry in this process and validating the results for silane thermal decomposition. Further extension of RMG to surface chemistry could allow for more rigorous comparison to CVD experiments.
■
° between high level calculations,39 Comparison of ΔfH298 CBS-QB3, G3//B3LYP, and group additivity. (ZIP)
REFERENCES
(1) Coltrin, M. E.; Kee, R. J.; Evans, G. H. A Mathematical Model of the Fluid Mechanics and Gas-Phase Chemistry in a Rotating Disk Chemical Vapor Deposition Reactor. J. Electrochem. Soc. 1989, 136. (2) Giunta, C. J.; McCurdy, R. J.; Chapple-Sokol, J. D.; Gordon, R. G. Gas-phase kinetics in the atmospheric pressure chemical vapor deposition of silicon from silane and disilane. J. Appl. Phys. 1990, 67, 1062−1075. (3) Ho, P.; Coltrin, M. E.; Breiland, W. G. Laser-induced fluorescence measurements and kinetic analysis of Si atom formation in a rotating disk chemical vapor deposition reactor. J. Phys. Chem. 1994, 98, 10138−10147. (4) Dollet, A.; De Persis, S.; Teyssandier, F. Rate constants of reactions involving SiH4 as an association product from quantum RiceRamsperger-Kassel calculations. Phys. Chem. Chem. Phys. 2004, 6, 1203−1212. (5) Dollet, A.; de Persis, S. Pressure-dependent rate coefficients of chemical reactions involving Si2H4 isomerization from QRRK calculations. J. Anal. Appl. Pyrolysis 2007, 80, 460−470. (6) Purnell, J. H.; Walsh, R. The Pyrolysis of Monosilane. Proc. R. Soc. London, Ser. A 1966, 293, 543−561. (7) Yuuki, A.; Matsui, Y.; Tachibana, K. A Numerical Study on Gaseous Reactions in Silane Pyrolysis. Jpn. J. Appl. Phys. 1987, 26, 747−754. (8) Frenklach, M.; Ting, L.; Wang, H.; Rabinowitz, M. J. Silicon particle formation in pyrolysis of silane and disilane. Isr. J. Chem. 1996, 36, 293−303. (9) Swihart, M. T.; Girshick, S. L. Thermochemistry and kinetics of silicon hydride cluster formation during thermal decomposition of silane. J. Phys. Chem. B 1999, 103, 64−76. (10) Wong, H. W.; Li, X.; Swihart, M. T.; Broadbelt, L. J. Detailed kinetic modeling of silicon nanoparticle formation chemistry via automated mechanism generation. J. Phys. Chem. A 2004, 108, 10122− 10132. (11) Emeleus, H.; Reid, C. 220. The pyrolysis of disilane and trisilane. J. Chem. Soc. 1939, 0, 1021−1030. (12) Loh, S. K.; Jasinski, J. M. Direct kinetic studies of SiH3 + SiH3, H, CCl4, SiD4, Si2H6, and C3H6 by tunable infrared diode laser spectroscopy. J. Chem. Phys. 1991, 95, 4914. (13) Takahashi, J.; Momose, T.; Shida, T. Thermal Rate Constants for SiH4 ⇌ SiH3 + H and CH4 ⇌CH3 + H by Canonical Variational Transition State Theory. Bull. Chem. Soc. Jpn. 1994, 67, 74−85. (14) Bhandarkar, U. V.; Swihart, M. T.; Girshick, S. L.; Kortshagen, U. R. Modelling of silicon hydride clustering in a low-pressure silane plasma. J. Phys. D: Appl. Phys. 2000, 33, 2731−2746. (15) Robertson, R.; Hils, D.; Chatham, H.; Gallagher, A. Radical species in argon-silane discharges. Appl. Phys. Lett. 1983, 43, 544−546. (16) Watanabe, Y. Contribution of short lifetime radicals to the growth of particles in SiH4 high frequency discharges and the effects of particles on deposited films. J. Vac. Sci. Technol., A 1996, 14, 995.
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.6b02402. libraries/: Chemkin formatted files and species dictionaries for rate parameters in reaction libraries. 2,5 NoPDep_base/: RMG generated mechanism for silane thermal decomposition without pressure dependence. PDep_base/: RMG generated mechanism for silane thermal decomposition including pressure dependence. PDep_no_radicals/: RMG generated mechanism for silane thermal decomposition with pressure dependence and generated without radical reaction families. PDep_y0/: RMG generated mechanism for silane thermal decomposition with pressure dependence valid over different initial silane concentrations. Residence Times.ipynb: An iPython notebook detailing calculations of error bars in Figure 5. spc_geom_energies/: Optimized geometries and energies for silicon hydride radicals. ts_geom_energies/: Optimized geometries and energies of transition states of silicon hydride hydrogen abstraction reactions. thermo_comparison_tables.csv: 12514
DOI: 10.1021/acs.iecr.6b02402 Ind. Eng. Chem. Res. 2016, 55, 12507−12515
Article
Industrial & Engineering Chemistry Research (17) Gao, C. W.; Allen, J. W.; Green, W. H.; West, R. H. Reaction Mechanism Generator: Automatic construction of chemical kinetic mechanisms. Comput. Phys. Commun. 2016, 203, 212−225. (18) RMG: Reaction Mechanism Generator. Open-source software, http://reactionmechanismgenerator.github.io. (19) Song, J.; Stephanopoulos, G.; Green, W. Valid parameter range analyses for chemical reaction kinetic models. Chem. Eng. Sci. 2002, 57, 4475−4491. (20) Benson, S. W.; Buss, J. H. Additivity Rules for the Estimation of Molecular Properties. Thermodynamic Properties. J. Chem. Phys. 1958, 29, 546. (21) Song, J. Building robust chemical reaction mechanisms: next generation of automatic model construction software. Ph.D. thesis. Massachusetts Inst. Technol. Dept. Chem. Eng.: Cambridge, MA, 2004. (22) Onischuk, A. A.; Strunin, V. P.; Ushakova, M. A.; Panfilov, V. N. Studying of silane thermal decomposition mechanism. Int. J. Chem. Kinet. 1998, 30, 99−110. (23) Adamczyk, A. J.; Reyniers, M.-F.; Marin, G. B.; Broadbelt, L. J. Exploring 1,2-hydrogen shift in silicon nanoparticles: reaction kinetics from quantum chemical calculations and derivation of transition state group additivity database. J. Phys. Chem. A 2009, 113, 10933−46. (24) Adamczyk, A. J.; Reyniers, M.-F.; Marin, G. B.; Broadbelt, L. J. Kinetic correlations for H2 addition and elimination reaction mechanisms during silicon hydride pyrolysis. Phys. Chem. Chem. Phys. 2010, 12, 12676−96. (25) Adamczyk, A. J.; Reyniers, M.-F.; Marin, G. B.; Broadbelt, L. J. Kinetics of substituted silylene addition and elimination in silicon nanocluster growth captured by group additivity. ChemPhysChem 2010, 11, 1978−94. (26) Baboul, A. G.; Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-3 theory using density functional geometries and zero-point energies. J. Chem. Phys. 1999, 110, 7650. (27) Wu, S. Y.; Raghunath, P.; Wu, J. S.; Lin, M. C. Ab initio chemical kinetic study for reactions of H atoms with SiH4 and Si2H6: comparison of theory and experiment. J. Phys. Chem. A 2010, 114, 633−639. (28) Manion, J. A.; Huie, R. E.; Levin, R. D.; Burgess, D. R., Jr.; Orkin, V. L.; Tsang, W.; McGivern, W. S.; Hudgens, J. W.; Knyazev, V. D.; Atkinson, D. B.; Chai, E.; Tereza, A. M.; Lin, C.-Y.; Allison, T. C.; Mallard, W. G.; Westley, F.; Herron, J. T.; Hampson, R. F.; Frizzell, D. H. NIST Chemical Kinetics Database, NIST Standard Reference Database 17, Version 7.0 (Web Version), Release 1.6.8, Data version 2015.12. 2015; NIST: Gaithersburg, MD. (29) Zhao, Y.; Truhlar, D. G. A new local density functional for maingroup thermochemistry, transition metal bonding, thermochemical kinetics, and noncovalent interactions. J. Chem. Phys. 2006, 125, 194101. (30) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functions. Theor. Chem. Acc. 2008, 120, 215−241. (31) McLean, A. D.; Chandler, G. S. Contracted Gaussian basis sets for molecular calculations. I. Second row atoms, Z = 11−18. J. Chem. Phys. 1980, 72, 5639. (32) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R. Efficient diffuse function-augmented basis sets for anion calculations. III. The 3-21+G basis set for first-row elements, Li-F. J. Comput. Chem. 1983, 4, 294−301. (33) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-consistent molecular orbital methods 25. Supplementary functions for Gaussian basis sets. J. Chem. Phys. 1984, 80, 3265. (34) Curtiss, L. A.; Raghavach Ari, K.; Redfern, P. C.; Rassolov, V.; Pople, J. A. Gaussian-3 (G3) theory for molecules containing first and second-row atoms. J. Chem. Phys. 1998, 109, 7764−7776. (35) Zhao, Y.; Lynch, B. J.; Truhlar, D. G. Development and Assessment of a New Hybrid Density Functional Model for Thermochemical Kinetics. J. Phys. Chem. A 2004, 108, 2715−2719.
(36) Allen, J. W.; Green, W. H. CanTherm: Open-source software for thermodynamics and kinetics Included in: Reaction Mechanism Generator, v2.0.0. 2016; http://reactionmechanismgenerator.github.io. (37) Wong, H. W.; Nieto, J. C. A.; Swihart, M. T.; Broadbelt, L. J. Thermochemistry of silicon-hydrogen compounds generalized from quantum chemical calculations. J. Phys. Chem. A 2004, 108, 874−897. (38) Lay, T. H.; Bozzelli, J. W.; Dean, A. M.; Ritter, E. R. Hydrogen atom bond increments for calculation of thermodynamic properties of hydrocarbon radical species. J. Phys. Chem. 1995, 99, 14514−14527. (39) Katzer, G.; Ernst, M. C.; Sax, A. F.; Kalcher, J. Computational Thermochemistry of Medium-Sized Silicon Hydrides. J. Phys. Chem. A 1997, 101, 3942−3958. (40) Allen, J. W.; Goldsmith, C. F.; Green, W. H. Automatic estimation of pressure-dependent rate coefficients. Phys. Chem. Chem. Phys. 2012, 14, 1131−55. (41) Goodwin, D. G.; Moffat, H. K.; Speth, R. L. Cantera: An Objectoriented Software Toolkit for Chemical Kinetics, Thermodynamics, and Transport Processes. 2016; http://www.cantera.org. (42) Crosby, L. D.; Kurtz, H. A. Application of Electronic Structure and Transition State Theory: Reaction of Hydrogen With Silicon Radicals. Int. J. Quantum Chem. 2006, 106, 3149−3159. (43) Wang, W.; Feng, S.; Zhao, Y. Quantum instanton evaluation of the thermal rate constants and kinetic isotope effects for SiH4 + H → SiH3 + H2 reaction in full Cartesian space. J. Chem. Phys. 2007, 126, 114307. (44) Yu, X.; Li, S.-M.; Li, Z.-S.; Sun, C.-C. Direct Ab Initio Dynamics Studies of the Reaction Paths and Rate Constants of Hydrogen Atom with Germane and Silane. J. Phys. Chem. A 2000, 104, 9207−9212. (45) Arthur, N. L.; Miles, L. A. Rate constants for H + (CH3)4−nSiHn, n = 1−4. Chem. Phys. Lett. 1998, 282, 192−196. (46) Arthur, N. L.; Potzinger, P.; Reimann, B.; Steenbergen, H. P. Reaction of H-Atoms with Some Silanes and Disilanes - Rate Constants and Arrhenius Parameters. J. Chem. Soc., Faraday Trans. 2 1989, 85, 1447−1463. (47) Jasinski, J. M.; Chu, J. O. Absolute rate constants for the reaction of silylene with hydrogen, silane, and disilane. J. Chem. Phys. 1988, 88, 1678.
12515
DOI: 10.1021/acs.iecr.6b02402 Ind. Eng. Chem. Res. 2016, 55, 12507−12515
