Exact quantum dynamics study of the H(2S)+SiH+(X1Σ+)reaction on a new potential energy surface of SiH+2(X2A1)∗

Wen-Li Zhao(赵文丽) Rui-Shan Tan(谭瑞山) Xue-Cheng Cao(曹学成)Feng Gao(高峰) and Qing-Tian Meng(孟庆田)

1School of Information Science and Engineering,Shandong Agricultural University,Taian 271018,China

2Department of Science,Shandong Jianzhu University,Jinan 250101,China

3School of Physics and Electronics,Shandong Normal University,Jinan 250358,China

Keywords: Chebyshev wave packet method,reaction probability,integral cross section,rate constant

1. Introduction

During the past decades, the construction of accurate global potential energy surface (PES) of the methylene ion CH+2and related dynamics investigation have aroused great interest because of their importance in atmospheric chemistry and astrophysics.[1–6]Comparatively speaking,the global PES of SiH+2system and related dynamical calculations are seldom reported. For the silicon semiconductor industry, SiH+and H are known to be two reactive intermediate species in chemical-vapor deposition process, and the investigation on H+SiH+reaction is of great practical significance.[7]Furthermore,SiH+and H are known to be abundant in interstellar space, and the study of H+SiH+reaction can help us understand the physical and chemical processes taking place in astrophysics.[8]In 2019, our research team[9]firstly reported an accurate global ground state PES(hereafter GMS-PES)for SiH+2(X2A1)at multireference configuration interaction level including Davidson corrections by the many-body expansion method.Using the coupled states(CS)approximation,we also carried out Chebyshev wave packet dynamical calculations for the reaction

However,the CS approximation has been proved to be inaccurate in some chemical reactions due to neglecting the Coriolis coupling (CC) effect.[10,11]In order to obtain accurate dynamical results and quantitatively investigate the influence of CC effect on the H+SiH+reaction, an accurate Chebyshev wave packet calculation considering full CC effect is performed based on GMS-PES of SiH+2(X2A1)in this work.

From the theoretical point of view, with multireference configuration interaction (MRCI) calculation, a three-dimensional potential energy function for ˜X2A1and ˜A2B1states of SiH+2was constructed in 1994 by Baueret al.[14]They also calculated the transition wavenumbers between ˜X2A1and ˜A2B1states based on the potential energy function, and the theoretical and experimental molecular parameters agree well. In 2002, Vachet al.[15]performed a molecular dynamics simulation to obtain the reactive cross sections of Si+(2P)+H2reaction in the collision range of 0.3 eV to 10 eV based on a semi-empirical PES of SiH+2system. In 2004, Chaˆabaneet al.[16]carried out an investigation of molecular dynamics simulation with semi-empirical quantum chemistry on Si+(2P)+H2reaction in the energy range of 0.5 eV and 10 eV,and the reaction cross sections for different rotational and vibrational excitation of H2(v0,j0)and reaction mechanisms were discussed,which indicated the final energy distributions of the SiH+products were controlled by the initial translational,rotational and vibrational energies of H2.

In this paper, based on the GMS-PES of SiH+2(X2A1)system,[9]we carry out an accurate quantum dynamical calculation of H+SiH+reaction by using the Chebyshev wave packet method with considering full CC effect. Chebyshev wave packet can be used to study the systems with large mass due to the more advanced scaling laws over time-independent methods.[17–21]As the wave packet propagates in real space without approximation and reaction properties of all energy points can be calculated by a single wave packet propagation,the Chebyshev wave packet propagation scheme used in the present paper have been proved to be accurate and efficient.[18]The paper is organized as follows: the theoretical and computational methods used in this paper are described in Section 2.Section 3 presents the main results of the accurate quantum dynamics with considering full CC effect including reaction probabilities, integral cross sections (ICS) and rate constants for the title reaction and discusses in detail the influence of CC effect on the SiH++H reaction. Finally,Section 4 concludes the present work.

2. Theoretical approach

In the reactant(H+SiH+)Jacobi coordinatesR(distance of H from the center of mass of SiH+),r(internuclear distance of SiH+)andγ(the angle betweenRandr),the Hamiltonian of SiH+2(X2A1)system can be described as

where the quantityµr=mSimH/(mSi+mH) is the reduced mass of the SiH+,µR=mH(mSi+mH)/(2mH+mSi) is the three-body reduced mass of the H+SiH+system, andV(R,r,γ) is the interacting potential energy of SiH+2(X2A1)system which can be obtained from PES. The rotational angular momentum operator and the orbital angular momentum operator of SiH+are denoted as ︿jand ︿L, respectively, which can be expressed as

In Eq.(3), ˆJdenotes the total angular momentum operator, ˆJzand ˆjzare the body-fixedzcomponents of total angular momentum operator and diatomic angular momentum operator,and ˆJ+(ˆJ−) and ˆj+(ˆj−) are the corresponding raising (lowering) operators. The last two terms in Eq. (3) are CC effect terms,which are ignored in CS approximation.

As neither discretization nor interpolation errors are introduced,the stable and accurate recursion relationship in the Chebyshev wave packet scheme are employed[18–21]

where ˆHnorm=( ˆH −H+)/H−is the normalized Hamiltonian withH±=(Hmax±Hmin)/2,in whichHmax(Hmin)is the upper (lower) bounds of spectral range. The wave packets are multiplied with Gaussian shaped damping function at the grid edges to enforce the outgoing boundary condition,

In Eq.(6),Nis normalization constant,k0,R0andδare initial momentum,central position and width of initial wave packet,respectively,and|ϕi〉is a specific rovibrational eigenfunction.

By using a flux-based approach,[22]the initial statespecified total reaction probability is expressed as

whereθ ≡arccos((E −H+)/H−)is the Chebyshev angle,rfdefines the dividing surface in the product channel andai(E)is the amplitude in the initial wave packet, which can be expressed as

The initial state-specified ICS is calculated by summing up the reaction probabilities for all quantum numbersJ’s of the total angular momentum for a given collision energyEc

wheref,kB,andTare the degeneracy factor,Boltzmann constant and the temperature,respectively.

3. Result and discussion

3.1. Potential energy surface and minimum energy path

Fig. 1. Two contour plots for different configurations. (a) Contour lines for bending angle ∠[H–Si–H]+=90◦figuration starting from −5.82 eV.(b)Contour lines for bending angle ∠[H–Si–H]+=180◦starting from −4.63 eV.Both tow panels are drawn with an increment of 0.20 eV.

Fig.2. The MEP of the different approaching angles(∠[H–Si–H]+=90◦&180◦)for H+SiH+→Si++H2 reaction.

3.2. Reaction probabilities

By using the Chebyshev wave packet method, the accurate quantum dynamical investigation is employed on H(2S)+SiH+(X1Σ+)→Si+(2P)+H2(X1Σ+g) reaction based on the GMS-PES for SiH+2(X2A1) system. Through a systematic test of the convergence withJ= 0, the optimum computational parameters in Jacobi coordinates are determined,which can be found in Ref.[9]. Due to the existence of deep potential well in reaction path and huge mass of SiH+2(X2A1)system, the large computational parameters and plenty of partial waves bring out a huge challenge to our quantum dynamical calculations,[9]especially for the exact quantum calculations considering full CC effect. In order to get the exact ICSs and rate constants of the title reaction, the probabilities of all the partial waves up toJ=90 are calculated in the collision energy of 0 eV–1.0 eV. The number of partial waves in CS approximation is 100,[9]which is larger than the corresponding ones of exact calculations in the same collision energy range. As seen in Eqs. (2) and (3), the increment of centrifugal potential barrier becomes faster with the increase ofJwhen the CC effect term is introduced into the Hamiltonian,which will reduce partial waves for exact quantum calculation in the same collision energy. In this work,the expensive computations are only carried out on probabilities ofJ=0,5,10,...,90 for the title reaction,and the others are obtained by cubic spline interpolation method.[23]

The reaction probabilities ofJ=0,20,40,60 and 80 with exact quantum calculations are displayed in Fig.3 for collision energy up to 1.0 eV. The distinctive feature is zero threshold for the probability ofJ= 0, which is owing to the barrierless insertion pathway of exothermic reaction. The reaction threshold shifts towards higher energy asJincreases,which is consistent withJ-shifting method[24]or capture model.[25,26]As so far, there are no experimental and other theoretical results of H+SiH+reaction compared with the results of this work, however, H+CH+reaction has been investigated extensively.[2,3,5,27,28]Furthermore, the topographical properties of CH+2(12A′) are essentially same with SiH+2(X2A1)system,[4,9]hence,the variation tendency and resonance structure of reaction probability,ICS and rate constant of H+SiH+reaction are compared with the ones of H+CH+reaction.

As the same with the H+CH+reaction, the strong oscillation structures of reaction probabilities are displayed in Fig. 3, which can be attributed to the formation of a SiH+2complex supported by stationary points of the SiH+2(X2A1)PES.[9]The reaction probability oscillates strongly in the low energy range and the amplitude of oscillation is reduced with increment of collision energy presumably as the lifetimes of resonance states become shorter. Compared with H+CH+reaction, the oscillation amplitude of reaction probabilities for H+SiH+reaction is relatively smaller,[5]which is because that the potential well depth of SiH+2system is shallower than that of CH+2and the resonance amplitude is mainly determined by the potential well depth.[5,9]For example, the maximum depth of CH+2system in C2Vsymmetry is 4.74 eV,[9]which is larger than the potential well depth of SiH+2system mentioned above. Although these two reactions are exoergic, the energy released for H+SiH+reaction is 1.32 eV,which is bigger than the exoergicity(0.496 eV)from the reactant H+CH+to the product C++H2. Therefore, the reaction probabilities of H+SiH+reaction are larger than the ones of H+CH+reaction overall.

As displayed in Fig.4, the reaction probabilities forJ=10, 20, 30, 40, 50, and 60 of the exact quantum calculations with considering full CC effect are compared with the corresponding results of CS approximation. Obviously, both CC and CS results show sharp oscillating structures in the collision energy of 0 eV–1.0 eV for H+SiH+reaction due to the deep well of SiH+2system. However, the oscillation amplitudes of reaction probabilities for CC calculations are much smaller than the corresponding CS results.This is because that the CC effect dominated by long-range forces will increase the collision time by orientating the reactants and decrease the lifetime of resonance states.[11,29]For lowJ’s(J=10 and 20),the CC probabilities are larger than the CS ones at low collision energy(aboutEc≤0.75 eV)and become slightly smaller than the CS results with the increment of collision energy.For highJ’s (J= 30, 40, 50 and 60), the CC probabilities are smaller than the CS ones at low collision energy. WhenJ=30, 40 and 50, the CC results are slightly larger than the CS ones at high collision energy(aboutEc≥0.70 eV).WhileJ=60,the CC results are much larger than the CS ones above 0.55 eV. Overall, the differences between the CS and CC results are observed to be remarkable in the region where the reaction probabilities of CS fluctuate greatly. While,in the region with small oscillation of CS probabilities, the results of CC are basically consistent with those of CS. In CS approximation,SiH+molecule is confined tok=0(the projection of total angular momentumJonto the body-fixedzaxis) quantum state,which will decrease the collision time,increase the resonances and lead to overestimation or underestimation of the reaction probability.[11,29,30]

Fig.3. The exact reaction probabilities varied with collision energy of 0 eV–1.0 eV for J=0,20,40,60 and 80.

Fig.4. The reaction probabilities of CC and CS calculations for H+SiH+→Si++H2 reaction at J=10,20,30,40,50,and 60.

3.3. Integral cross section and rate constant

The ICSs of CC and CS calculations are plotted over the collision energy range of 0 eV–1.0 eV in Fig. 5. It is obvious that both ICSs have no threshold. They decrease sharply with the increase of collision energy in the low energy region and then decrease slowly as collision energy increases. The trend of ICS for H+SiH+reaction varied with collision energy is similar to H+CH+reaction, which is determined by reaction mechanism of an exothermic reaction without potential barrier. After a weighted sum of all reaction probabilities,the oscillatory structure of ICS for CS approximation is still relatively obvious. While for CC ICS, the oscillatory structure has been smeared out mostly. As discussed in Ref.[29],CC effect has a significant influence on resonance by breaking the resonance complex. In the exact quantum calculation with considering full CC effect,ICS becomes smoother and less oscillating compared with the CS ones. In the collision energy range of 0 eV–0.05 eV, the ICSs of CC are much larger than CS ones. When 0.05 eV≤Ec<0.1 eV,CS and CC ICSs are roughly equal. When 0.1 eV≤Ec<0.5 eV,the CC results are slightly smaller than CS ones.When 0.5 eV≤Ec<0.7 eV,the CC results are basically same with CS ones. When 0.7 eV≤Ec≤1.0 eV,the CC results are slightly larger than CS ones.

According to Eq. (10), the initial state specific (vi=0,ji= 0) thermal rate constants for H(2S)+SiH+(X1Σ+)→Si+(2P)+H2(X1Σ+g) reaction are calculated from the energy resolved ICSs. The rate constants varied with temperature range from 5 K to 1000 K for both CC and CS calculations are displayed in Fig. 6. The rate constant of CC calculations increases sharply in the temperature range from 5 K to 75 K and then creeps up in the range of 75 K–175 K. From 175 K up to 1000 K,the CC result decreases slightly and are nearly constant. While for CS rate constant, a monotonic increasing behavior is shown. Overall, the CS result increases rapidly in the temperature range from 5 K to 100 K and then grows gradually above about 100 K.By comparison, CC and CS results are well consistent below 16 K. In the range of 16 K–100 K, CC rate constant increases more rapidly than CS ones, and the deviation of them becomes larger with the temperature increasing and reaches the maximum value about at 100 K. For instance, whenT= 25 K, the CC rate constant isk=9.2×10−10cm3·molecule−1·s−1the CS result isk=7.0×1010cm3·molecule−1·s−1, and the relative error is 24.1%. WhenT=100 K, the CC and CS rate constants are 1.7×10−9and 9.4×1010cm3·molecule−1·s−1, respectively,and the relative error is 44.4%. In the range of 100 K–1000 K,the CC rate constant decreases slightly while the CS ones gradually increase, and the relative error of them obviously decreases. Compared with the reaction of H+CH+,the rate constants of these two reactions are in the same order of magnitude at the same temperature, and further observation shows that the varied trend of CC rate constants of H+SiH+reaction is more similar to the ones of H+CH+collision.[3,5]

From the ICSs and rate constants discussed above,it can be seen that the deviation of CC and CS ICSs is small in the most of collision energy range except for the range 0 eV–0.05 eV. But for rate constants, the deviation of them is considerable in the temperature range of 16 K–1000 K. In order to quantitatively investigate the effect of ICSs in the range of 0 eV–0.05 eV on the rate constants in the range of 0 K–1000 K,according to Eq.(10)the function of proportionPr(T)varied with temperatureTis calculated by

where the numerator is the integration of partition functionEce−Ec/kBTover collision energyEcbelow 0.05 eV,while the denominator is the total integration of 0.00001 eV–1 eV.In order to ensure the convergence of integral calculation,the lower limit of integration is 0.000001 instead of zero.The calculated function of proportionPr(T)varied with temperatureTis displayed in Fig.7.

Fig.5. The ICSs of CS and CC calculations for H+SiH+→Si++H2 reaction over collision energy of 0 eV–1.0 eV.

Fig.6. The initial state(ν=0, j=0)specified rate constant of H+SiH+→Si++H2 reaction for CC and CS calculations obtained on the GMS PES.

In the above figures, the effect of ICSs in the range of 0 eV–0.05 eV on the rate constants in the range of 0 K–1000 K can be described as follows: when 0

Fig.7. The function of proportion Pr(T)varied with temperature T.

4. Conclusion

Based on the global accurate GMS-PES of SiH+2(X2A1),the Chebyshev quantum wave packet method is employed to study the dynamics of H(2S)+SiH+(X1Σ+)→Si+(2P)+H2(X1Σ+g) reaction in this work, and the initial state specified (ν=0,j=0) probabilities, ICSs and rate constants are provided. The dynamical behaviors of these calculated results are according with an exothermic reaction without potential barrier. In order to investigate the effect of CC on the title reaction,a large number of partial waves(J ≤90)in the collision energy range of 0 eV–1.0 eV are included, which brings out a huge challenge to exact quantum calculations with considering full CC effect. By comparing the CC probabilities with the corresponding CS ones,it is found that the CC effect dominated by long-range forces will increase the collision time by orientating the reactants and decrease the lifetime of resonance states,which means that neglecting CC effect will decrease the collision time,increase the amplitude of oscillation and lead to overestimation or underestimation of the reaction probability.For ICSs and rate constants, it is found that the deviation of CC and CS ICSs is small in the most of collision energy range except for the range of 0 eV–0.05 eV, while the deviation of both rate constants is considerable in the temperature range of 16 K–1000 K, which can be quantitatively illustrated by the function of proportionPr(T). In conclusion,the CC effect is of great importance and can not be ignored for the title reaction.