Silane-initiated nucleation in chemically active plasmas: validation of density functionals, mechanisms, and pressure-dependent variational transition state calculations

The growth of anionic silicon hydride clusters is a critically important process in nanodusty plasmas. In the current study, we focus on the formation of homologs of silylene (Sin+1H2n+2 , n = 3, 4) and silyl (SinH2n+1 , n = 4, 5) anions via anion–neutral reaction pathways. Species like silyl or silylene anions and their related elementary reactions, which are involved in the formation of silicon hydride clusters, were not used in developing exchange–correlation (xc) density functionals (i.e., they were not included in the training set of semiempirical density functionals); therefore, we explored the accuracy of various widely used xc density functionals based on reaction energies and barrier heights. Among the 21 density functionals we tested, M06-2X has the best performance for a hybrid functional, and MN15-L has the best performance for a local functional. Thermal rate constants of the elementary reactions involved in the reaction mechanism are calculated using M06-2X and multistructural canonical variational transition state theory with the small-curvature tunneling approximation (MS-CVT/SCT). The pressure dependence of unimolecular isomerization reactions is treated with systemspecific quantum RRK theory (SS-QRRK) and the Lindemann–Hinshelwood mechanism.


Introduction
The growth of nanoparticles in nanodusty plasmas is an active research field in plasma physics, chemistry, and engineering. Many physical and chemical processes are involved in the formation of nanoparticles in chemically active plasmas, including nucleation, isomerization, electron capture and ionization, and mass and momentum transport. Anion-neutral reactions are one of the major chemical reactions in silicon hydride clustering in silanecontaining reactive plasmas. 1,2 Rate constants of anion-neutral reactions are used in building transport equations or boundary conditions in many simulation studies for investigating the distribution of the sizes of nanoparticles and their population in a plasma. [3][4][5][6][7] However, due to the complexity of the systems and the unusual conditions of the chemical reactions (i.e., plasma), accurate experimental measurements for these anion-neutral reactions are very difficult and hence are scarce. The majority of the reaction rates needed for modeling nanodusty plasmas are empirically estimated, and Agarwal and Girshick 8 noted that ''there is considerable uncertainty in the values of rate constants for anion-neutral reactions that are primarily responsible for clustering in this system, and thus any correction factor one might apply for the predicted nucleation rate would itself be highly uncertain''. Theoretical calculations of chemical structures, energetics, and thermal rate constants can play an important role in reducing the uncertainty involved in plasma modeling. Although fast empirical methods for estimating thermodynamic functions and rate constants exist, 9 we do not use such approach in the current work; thermodynamic functions and thermal rate constants reported in this work are computed from ab initio calculations.
In the present study, we focus on clusters that contain no more than five silicon atoms. The chemical mechanism of the initial step of the clustering process, i.e., the formation of (H 3 Si) 2 HSi: À and H 3 SiH 2 SiHSi À , has been investigated in a previous work. 10 In the current work, since the studied kind of system is not represented in the training set of any density functional known to us, we carry out benchmark tests for various density functionals, and this provides an opportunity to investigate the transferability of various semiempirical exchange-correlation functionals. Chemical structures, energetics, and thermal rate constants in the mechanisms are computed using the best functional among the functionals we tested in this work.
Two unimolecular isomerization reactions are involved in the reaction mechanism we propose in this work, and the pressure dependences of their thermal rate constants are estimated using system-specific quantum Rice-Ramsperger-Kassel (SS-QRRK) theory combined with the Lindemann-Hinshelwood thermal activation mechanism. The SS-QRRK method was proposed recently 11 as a convenient way to use variational transition state theory to treat pressure dependences of chemical reaction rates, and it was applied to a chemical activation mechanism; the present article is the first application to a thermally activated unimolecular reaction. For a chemically activated unimolecular reaction, the low-pressure rate constants of the formation of the stabilized adduct are lower than the high-pressure-limit, and the rate constants of the further isomerization/dissociation of the formed adduct are larger than the high-pressure equilibrium rate constant; and in thermally activated unimolecular reactions, the low-pressure rate constants are smaller than the high-pressure equilibrium rate constant, and hence the pressure effect for this kind of reaction is called ''falloff.'' 2. Theoretical background 2.1. High-pressure-limit thermal rate constants High-pressure-limit thermal rate constants are computed using multi-structural canonical variational transition state theory with the small-curvature tunneling approximation (MS-CVT/SCT) as follows: [12][13][14][15] where F act is multi-structural torsional potential anharmonicity factor of activation 16,17 computed by MS-T(C) method; 18 k SCT 1 is a small-curvature tunneling transmission coefficient for the lowest-energy conformer of the transition state, G CVT 1 is CVT variational transmission coefficient for the lowest-energy path, and k TST 1 is the conventional transition state theory rate constant computed based on the lowest-energy path: where Q ‡ and Q R are rigid-rotor-harmonic-oscillator partition functions for the transition state structure and reactant; F R is the reactants partition function per unit volume. V ‡ is the barrier height, k B is the Boltzmann constant, h is Planck's constant, and T is absolute temperature. Expressing the TST rate constant in the form of eqn (2), the rotational partition function does not contain rotational symmetry number; the rotational symmetry number is included in the reaction symmetry number s rxn as follows: 19 where s R and s ‡ are the rotational symmetry number for the reactant and transition state structure, which are equal to the order of the rotational subgroup for polyatomic molecule; rotational symmetry number is 1 for heteronuclear diatomic molecule and 2 for homonuclear diatomic molecule. The reaction symmetry number we used here excludes the non-superimposable mirror-image conformers; the contributions from these mirrorimage conformers are treated within MS-T method. However, the contributions from enantiomers should be included in the reaction symmetry number or the rotational partition function if the MS-T method is not applied. 19,20 The reaction symmetry numbers for forward and reverse reactions involved in the current study are tabulated in the ESI. †

System-specific quantum RRK theory with Lindemann-Hinshelwood theory
In this section, we consider the effect of pressure on the unimolecular reactions by using SS-QRRK theory 11 with the Lindemann-Hinshelwood mechanism. For a unimolecular reaction, the following Lindemann-Hinshelwood thermal activation mechanism is considered: 21 AðTÞ þ M À ! À where A(T) represents the thermally equilibrated reactant at temperature T; M is the bath gas; A*(E) is the rovibrationally excited molecule with total rovibrational energy E; and P is the isomerization product. This mechanism includes the RRK assumption that energy in A* is rapidly statistically redistributed among modes subject only to the constraint of total energy E so that the reactivity of A* is simply a function of total energy. The rate constant of energization is k 1 , and the rate constant of isomerization is k 2 , and by the RRK assumption both of these rate constants depend on the total energy of A*; and k 1 also depends parametrically on temperature T. The rate constant of de-energization is k c and is treated as temperaturedependent but energy-independent, which is the strong collision assumption. However the strong-collision assumption is mitigated in the present work by computing k c as the product of the Lennard-Jones collision rate constant and a collision efficiency factor b c , where the latter is computed by using Troe's modified collision model. 22,23 The pressure-dependent unimolecular reaction rate constant k uni (T, p) for the above mechanism is: [24][25][26] where K(E,T) is the equilibrium constant of the first step, [M] is the concentration of bath gas, p is pressure, and E 0 is the threshold energy. Note that K(E,T) is mixed-ensemble equilibrium constant, representing the thermal equilibrium of species A* in a microcanonical ensemble at energy E with species A in a thermal ensemble with temperature T. Therefore, where r(E) is the density of states, and Q A (T) is the rovibrational partition function of A.
In QRRK theory, 24,27 one models the states of A and A* as the discrete states of a system with s uncoupled harmonic oscillator modes, each with frequency o in wave numbers (e.g., cm À1 ). Then the integrals over E are replaced by sums over n, where n is the number of quanta excited at energy E (n = E/hc o, where the zero of energy is the potential energy of the equilibrium structure of A). With this model and with [M] given by the ideal gas law, eqn (4) can be written as: where and R is ideal gas constant, and c is the speed of light. Notice that in the limit of p -N, the high-pressure-limit unimolecular rate constant is recovered in eqn (4) and (6), and in the lowpressure limit, the unimolecular rate ''constant'' is no longer a constant (i.e., no longer independent of concentrations), but rather is proportional to [M] or p.
In SS-QRRK, we use eqn (6)- (8) with o taken as the geometric mean of the s vibrational frequencies of A (where s is 3N À 6, and N is the number of atoms in A), and with E 0 and the microcanonical isomerization rate constant k 2 (E) parameterized by using QRRK theory calibrated to match the MS-CVT/SCT canonical rate constant at high pressure. In QRRK theory, the microcanonical rate constant is a frequency factor A (with units of reciprocal time) times the probability that a molecule with n quanta of vibrational theory has at least m quanta in one of the modes (the reactive mode), where m is given in terms of the threshold energy by eqn (8). This yields 24 Substituting this into eqn (6), carrying out the sum, and taking the limit of p -N yields the Arrhenius form: 19,24 Thus, to parameterize QRRK theory, the parameter m is calculated from the MS-CVT/SCT Arrhenius activation energy E MS-CVT/SCT a (T): and the frequency factor A in QRRK theory is set equal the MS-CVT/SCT Arrhenius pre-exponential factor, i.e., Note that both E 0 and A depend on temperature in the parameterized rate expression. Therefore k 2 (E) in eqn (6) becomes k 2 (E,T) given by eqn (9) and (11)- (13). The sum in eqn (6) is evaluated with a step size of one quantum, and the factorials therefore all have non-integer arguments and are evaluated using gamma functions. The de-energization rate constant is modeled as in our previous work, 11 using empirical Lennard-Jones parameters, the average energy transferred, and the energy dependence factor F E of the density of states (which is the thermal population of unimolecular states above the threshold energy of the reactant normalized by a density of states factor at the threshold energy); these quantities are given in Section 3.4.

Electronic structure calculations
Initial geometries of the reactants, products, and transition state structures are optimized with the M08-HX functional 28 and the MG3S basis 29 for benchmark studies; tight convergence criteria are used for both the SCF calculations and the geometry optimizations. For silicon atoms, the MG3S basis is equivalent to the 6-311+G(3d2f) basis, and for hydrogen atoms it is equivalent to the 311G(2p) basis. [30][31][32] All the density functional integrations are carried out with a grid of 99 radial shells around each atom and 974 angular points per shell. 33 All the electronic structure calculations are performed with a locally modified version of the Gaussian 09 software. 34,35 In order to test the accuracy of various popular density functionals for silyl or silylene anions reactions, a benchmark study was carried out based on classical barrier heights and classical energies of reactions (these are relative Born-Oppenheimer potential energies at stationary points and are exclusive of zero-point or thermal vibrational energy). Five reactions (as summarized in Table 1) were selected from our proposed reaction mechanisms for use in the benchmark study.
The reference energy values are computed by CCSD(T) 56 / CBS, where the complete basis set (CBS) limit is obtained by the following strategy: 57 where SBS means the small basis set (which is aug-cc-pVTZ 37,58 in the current work), and the MP2 59 /CBS energy is computed by respectively extrapolating the Hartree-Fock (HF) exchange energy and MP2 correlation energy as follows: 60-64 In the above equations, we use X = 4 (which is taken as aug-cc-pVQZ) and X À 1 = 3 (which is taken as aug-cc-pVTZ). Finally,

Direct dynamics calculations
For the direct dynamics calculations, all the species are re-optimized by M06-2X/MG3S, which was selected as the level for direct dynamics calculations based on the benchmark tests (as will be discussed in Section 4.1). Canonical variational transition state theory calculations with the small-curvature tunneling approximation were carried out in non-redundant internal coordinates 65,66 with a step size of 0.002 a 0 (note: 1 a 0 = 1 bohr = 0.5292 Å). The minimum energy paths (MEPs) are computed using Page-McIver algorithm 67 from À2.0 to +2.0 a 0 . The RODS algorithm 68 was used to re-orient the generalized-transition-state-theory dividing surface. A scaling factor 0.970 69 was used to scale all the vibrational frequencies in the generalized normal mode calculations.
For bimolecular reactions with a negative barrier, the smallcurvature tunneling transmission coefficient was computed using the ion-dipole complex as the initial state; the final bimolecular reaction rate constant is the product of the so-obtained tunneling transmission coefficient with the bimolecular reaction rate constant computed without tunneling. This is consistent with the fact that the tunneling calculation is performed for the high-pressure limit where the ion-dipole complex is fully thermalized (the issue of tunneling from the states of the precursor complex at lower energies than the bimolecular reactant ground state is discussed elsewhere 70,71 ). Multi-structural torsional anharmonicity (MS-T) rovibrational partition functions were computed using the MSTor program; 72 VTST calculations were performed with the Polyrate 73 and Gaussrate 74 programs.

Torsional anharmonicity
The multi-structural torsional anharmonicity (MS-T) rovibrational partition functions are computed based on coupled effective torsional potentials. The local periodicities of -SiH 3 groups are set to be 3. The MS-T partition functions include the contributions from all the distinguishable conformational structures including non-superimposable mirror images.

Parameters used in de-energization
In the current study, Ar gas, which is commonly used in chemical vapor decomposition 75 (CVD) and in studying nucleation processes in plasmas, 8 is selected as the bath gas used in estimation of the pressure-dependent rate constants of the unimolecular isomerization reactions. Lennard-Jones parameters e/k B and s are taken for Ar as 120 K and 3.4 Å 76 and for Si 4 H n species as 254 K and 5.8 Å, as used in previous silane plasma modeling. 77  Table 2 Exchange-correlation functionals tested in the current work and their percentage of non-local Hartree-Fock exchange (% X) The energy transfer parameter hDEi, which is the average energy transferred during both energization and de-energization processes and which is used for computing the collision efficiency factor, is chosen to be 740 cal mol À1 ; this value has been used previously for modeling SiH 4 colliding with Ar. 78 The Si clusters we considered in the current work are larger than SiH 4 , and therefore one might hypothesize that a larger energy transfer parameter hDEi should be used. To test this, we repeated the calculations for reaction RB, step 1 with a doubled energy transfer parameter of 1480 cal mol À1 , and we found that the obtained falloff curves (shown in Fig. S2 in ESI †) are not sensitive to this change. Doubling the energy transfer parameter leads to larger k(p) values (i.e., stronger collisions and smaller falloff effect) and at most a factor of 2 difference. (The maximum effect is at 1500 K, 0.001 bar.) The collision parameters and energy transfer parameters, in principle, could be determined from theoretical trajectory calculations; 79 we do not use such approach in the current work based on computational cost and the desire for a simple method that can be widely used in mechanism development.

Intramolecular hydrogen migration
The energy dependence factor F E of the density of states can be directly computed using Troe's definition: 22 with density of states computed from MS-T partition function by inverse Laplace transform; alternatively, F E can be computed by the empirical Whitten-Rabinovitch approximation using eqn (6) and (8) in Troe's work. 80 F E is used to compute the collision efficiency coefficient b c using the following equation: For falloff calculations on (SiH 3 ) 2 SiHSiH À , we used F E values computed from MS-T partition function in eqn (19)

Benchmark of various density functionals
Although coupled cluster theory with single and double excitations and quasi-perturbative connected triple excitations, i.e., CCSD(T), is often viewed as the gold standard in quantum chemistry, it uses a single configuration state function as the reference wave function (i.e., it is a ''single-reference'' method), and therefore it might not be appropriate to use CCSD(T) as a benchmark for systems with strong multi-reference characters (such as bond dissociation 81 and transition metal chemistry 82 ).
To ascertain its expected reliability, we computed the T 1 diagnostic 83 values for all the species involved in the reactions that are used for benchmark study, since this diagnostic has been proposed as a measure of the suitability of singereference coupled cluster theory. The molecule with the largest T 1 diagnostic value is the transition state structure of reaction (R1) (Si 4 H 10 À ), for which the value is 0.0252 computed by CCSD(T)/aug-cc-pVTZ; this value is much smaller than 0.045, which has been suggested 84-86 as a criterion for the applicability of single-reference methods to open-shell molecules, and hence we concluded that CCSD(T) can serve as a reference for testing other methods in this work. which is usually quoted as ''chemical accuracy''. M06-2X, M05-2X, M08-SO, and MPW1K are the four best performing functionals and all have overall MUEs within the chemical accuracy criterion. Among the three tested local density functionals (M11-L, MN12-L and MN15-L), MN15-L possesses the smallest overall MUE, which is 1.25 kcal mol À1 ; it is encouraging that this surpasses the accuracy of many popular hybrid density functionals (such as B3LYP and oB97X-D) and the recently developed MGGA_MS2h. This is encouraging for the very new MN15-L functional since there is no silyl or silylenerelated species in its training set and since local functionals are usually less accurate than hybrid functionals for reaction energies and barrier heights.
Note that the density functional calculations converge more rapidly with respect to basis set than CCSD(T), and so CBS extrapolations were not needed. A great advantage of density functional theory in this respect is that the reaction path calculations and the required Hessians along the reaction path are affordable even with the reasonably well-converged MG3S basis set.
Potential energy diagrams (relative Born-Oppenheimer potential energy E in kcal mol À1 with respect to reactants) for reaction mechanisms of reactions RA-RD are shown in Fig. 2-5. The elementary steps involved and their classical energies of reaction (DE) are listed in Table 3. Energetic values are computed based on geometries optimized at M06-2X/MG3S and single point energies calculated with the same method. All the elementary reactions are bimolecular reactions except for the first step of RB and the third step of RC, which are unimolecular isomerization reactions (intramolecular hydrogen transfer). We considered ion-dipole complex in four of the elementary bimolecular reactions, in which the barrier is negative: RA step 2, RB step 3, RC step 1, and RD step 1; the ion-dipole complex is used in computing tunneling transmission coefficient of the bimolecular reactions.

Multi-structural and torsional potential anharmonicity
Multi-structural and torsional potential anharmonicity factors of activation for the forward reactions are tabulated in Table 4, and for the reverse reactions are in ESI. † MS-T standard-state reaction enthalpy (DH ;MS-T rxn in kcal mol À1 ) and reaction Gibbs free energy (DG ;MS-T rxn in kcal mol À1 ) for all the reactions at 298 K and 1000 K are shown in Table 6, and their values at various temperatures are tabulated in ESI. † The numbers of distinguishable conformers found in the conformational searches for the species that have multiple structures are given in ESI. †  The errors for computing the thermal rate constants introduced by ignoring the multi-structural and torsional anharmonicity effects vary between different reactions and temperatures. For reactions RA step 1, RB step 3 and RC step 2, the errors are only slightly temperature dependent; the averaged errors (over 298-1500 K temperature range) for thermal rate constants computed without including MS-T effects for these reactions are respectively 32%, 43%, and 44%, with a standard deviation of 4%, 3%, and 5% between various temperatures. For reaction RB step 1, the MS-T effect can be ignored for temperatures below 800 K, at which temperatures the errors are all smaller than 15%; at 1500 K, one would have an error of 42% in thermal rate constant if the MS-T factor were not included. For reaction RC step 1 at 1500 K, the multistructural CVT/SCT rate constant is a factor of 5 smaller than the single-structural CVT/SCT rate constant; the ratio of the single-structural harmonic oscillator rovibrational partition function (SS-HO) to the MS-T rovibrational partition function is 0.138 for the TS and is 0.696 for the reactant (SiH 3 ) 2 SiH À , while the ratio of MS-HO partition function to SS-HO partition function is 2 for the TS and is 1 for (SiH 3 ) 2 SiH À , which indicates stronger couplings between torsional modes in the TS than in (SiH 3 ) 2 SiH À .

High-pressure-limit thermal rate constants
Calculated MS-CVT/SCT rate constants for all the forward reactions are listed in Table 5, and for reverse reactions in the ESI. † The following equations are used to fit the thermal rate constants: 70,87 k ¼ where A, n, E and T 0 are fitting parameters, T is temperature in Kelvin, and R is the ideal gas constant (1.9872 Â 10 À3 kcal mol À1 K À1 ). Fitting parameters for forward reactions are shown in Table 6, and those for the reverse reactions are in the ESI. † Fig. 6 shows the computed small-curvature tunneling (SCT) transmission coefficients at various temperatures for the third step in reaction RC. At low temperature, the thermal rate constant can be increased by an order of magnitude due to Fig. 4 Potential energy diagram for reaction mechanism of reaction RC; relative energies (E in kcal mol À1 , with respect to reactants (SiH 3 ) 2 SiH À + SiH 4 ) are computed at M06-2X/MG3S level.

Fig. 5
Potential energy diagram for reaction mechanism of reaction RD; relative energies (E in kcal mol À1 , with respect to reactants (SiH 3 ) 3 Si À + SiH 4 ) are computed at M06-2X/MG3S level. Table 3 Elementary steps in the proposed reaction mechanisms for reactions RA-RD, and their classical energies of reactions (kcal mol À1 ) at M06-2X/MG3S level

Reaction
Step quantum tunneling. The SCT tunneling transmission coefficient decays rapidly as temperature increases; k SCT decreases from 11.1 at 298 K to 2.0 at 500 K, and to 1.1 at 1500 K.

High-pressure-limit activation energy
The Arrhenius activation energies are defined by E a ¼ ÀR dðln kÞ dð1=TÞ (22) and are obtained by putting eqn (21) into eqn (22), which yields Computed activation energies for both forward and reverse reactions are shown in ESI. † The temperature dependence of activation energies of forward reactions is depicted in Fig. 7. Step 1 Step 2 Step 1 Step 2 Step 3 Step 1 Step 2 Step 3 Step 1 Step 2 298 9.08 Â 10 À32 3.71 Â 10 À10  Step 1 Step 2 Step 1 Step 2 Step 3 Step 1 Step 2 Step 3 Step 1 Step  Table 6 Fitting parameters for MS-CVT/SCT rate constants for forward reactions computed by M06-2X/MG3S and MS-T reaction enthalpy and reaction Gibbs free energy at 298 K and 1000 K a Step RA RB RC RD Step 1 Step 2 Step 1 Step 2 Step 3 Step 1 Step 2 Step 3 Step 1 Step 2 The temperature dependence of activation energy computed based on MS-CVT/SCT rate constants could be quite different from the one computed from conventional transition state theory (TST), because of the effects of tunneling, recrossing, and multi-structural torsional anharmonicity. For instance, in step 2 of reaction RC, the MS-CVT/SCT activation energy is 7.4 kcal mol À1 at 298 K, and it increases by 12.5 kcal mol À1 from 298 K to 1500 K; the conventional single-structural TST (without tunneling) activation energy at 298 K is 10.9 kcal mol À1 , which is 3.5 kcal mol À1 higher than MS-CVT/SCT activation energy at 298 K. We also plot the SS-TST, SS-TST/SCT, MS-CVT, and MS-CVT/SCT rate constant for RC step 2 at various temperatures in Fig. 8. The SS-TST/SCT natural logarithm rate constant curve decreases slower than the straight-line-shape of SS-TST curve as 1000 K/T increases; at high temperature, tunneling is negligible and SS-TST overlaps with SS-TST/SCT curve, while at low temperature, SS-TST curve is significantly lower than SS-TST/SCT curve. The MS-CVT/SCT curve differs negligibly from the SS-CVT/SCT curve due to the small MS-T effect.

Falloff effects for unimolecular isomerization reactions
Step 1 of reaction RB and step 3 of reaction RC are unimolecular isomerization reactions. The predicted falloff curves for these two reactions at various temperatures (in K) and pressures    (in bar) are shown in Fig. 9. Falloff curves are plotted as log 10 [k(p)/k(p = N)] versus pressure, where k(p) is the thermal unimolecular rate constant at pressure p, and k(p = N) is the high-pressure-limit rate constant computed by MS-CVT/SCT theory. The falloff curves for pressures from 10 À5 to 1 bar are shown in ESI. † At low and middle temperatures (T o 600 K), falloff effects for these two unimolecular reactions are negligible. At 600 K, the rate constant of RB step 1 at 1000 bar is 2.83 Â 10 2 s À1 , and at 0.01 bar it is 2.66 Â 10 2 s À1 ; the rate constant of RC step 3 at 600 K and 1000 bar is 2.19 Â 10 3 s À1 , while at the same temperature at 0.01 bar it is 1.90 Â 10 3 s À1 . For RB step 1, the collision efficiency coefficient b c at 298 K, 600 K, and 1000 K are 0.34, 0.14 and 0.03, respectively; and for RC step 3, collision efficiency coefficients at 298 K, 600 K, and 1000 K are 0.33, 0.12 and 0.02, respectively.
At high temperature, falloff effects become more significant. At 1500 K, the rate constant of RB step 1 is 1.45 Â 10 8 s À1 at 1000 bar and is 2.66 Â 10 6 s À1 at 1.0 bar, so the 1 bar result is a factor of 0.014 smaller than the high-pressure-limit. The rate constant of RC step 3 at 1500 K is 2.82 Â 10 8 s À1 at 1000 bar and is 1.25 Â 10 6 s À1 at 1 bar, so the high-pressure-limit rate constant is 574 times larger than the one at 1.0 bar.

Pressure-dependent activation energy
The activation energies for the unimolecular isomerization reactions RB step 1 and RC step 3 depend on both temperature and pressure. We show the pressure-dependent activation energies at 600 K, 800 K, and 1000 K in Fig. 10. At the highest pressures shown (10 3 bar), the activation energies are at the high-pressure limit. At 600 K and around 0.1 bar, where falloff effects on the rate become notable, the activation energies also start falling as the pressures decreases. At 800 K and 1000 K, the activation energy is significantly pressure dependent; the activation energy decreases almost linearly with respect to the pressure. At the transition pressure p 1/2 , 25 which is defined as the pressure at which the unimolecular rate constant is half of the high-pressure-limit value, the ratio of E a ( p = p 1/2 ) to E a ( p = N) is 0.6 for RB step 1 ( p 1/2 = 0.01 bar) and 0.6 for RC step 3 (p 1/2 = 0.03 bar) at 800 K, and it is 0.5 for RB step 1 ( p 1/2 = 0.3 bar) and 0.4 for RC step 3 ( p 1/2 = 1 bar) at 1000 K.

Summary
In the current work, we tested various exchange-correlation density functionals for an important system in nanodusty plasmas, in particular silylene and silyl anions reacting with silane molecules. Among the functionals we tested in this work, M06-2X is the most successful hybrid functional and MN15-L is the most successful local functional. Reaction mechanisms for the growth of silicon hydride clusters have been proposed, and the thermal rate constants of the elementary reactions involved in the reaction mechanisms were computed using multi-structural canonical variational transition state theory with the small-curvature tunneling approximation. Two unimolecular isomerization reactions are involved in the reaction mechanism, and their pressure dependent thermal rate constants were estimated based on system-specific quantum RRK (SS-QRRK) theory combined with Lindemann-Hinshelwood theory.
This work provides guidance for choosing density functionals for studying anion-neutral reactions in the silane-based reactive plasma. The methodology for computing thermal rate constants, whose values are rarely available experimentally, is also useful in estimating input kinetic data in plasma modeling for engineering applications. The extension of the SS-QRRK method to unimolecular isomerizations is also of more general use; it may be applied, for example, to atmospheric chemistry and combustion chemistry. Fig. 10 Predicted pressure-dependent activation energies (kcal mol À1 ) for unimolecular isomerization reactions RB step 1 and RC step 3 at 600 K, 800 K, and 1000 K. Solid lines are for reaction RB step 1; dashed lines are for reaction RC step 3.