Reaction mechanism of D+ND→N+D2 and its state-to-state quantum dynamics?

Ting Xu(許婷),Juan Zhao(趙娟),Xian-Long Wang(王憲龍),and Qing-Tian Meng(孟慶田),?

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

2 College of Science,Shandong Jiaotong University,Jinan 250357,China

Keywords:state-to-state quantum dynamics,time-dependent wave packet,D+ND,differential cross section

1.Introduction

The nitrogen(N)atom is one of the most abundant elements in the atmosphere.Furthermore,the reactions containing N atom play important roles in atmospheric chemistry,combustion chemistry and some explosion processes.[1-3]Consequently,the reactions of the N atom especially with the H atom and their isotopes,have attracted much attention from both theoristsand experimentalists.[4-19]Experimentally,Koshi et al.[4]and Davison and Hanson[5]directly detected N atoms by the atomic resonance absorption technique,and studied the N+H2reaction at the high temperatures(1950 K-2850 K)in a shock tube.Furthermore,the rate coefficients of the reactions NH+H[6]and NH+D[7]were given in a quasistatic laser- flash photolysis,laser-induced fluorescence system at room-temperature.

The theoretical investigation of the related reaction dynamics,such as with time-dependent wave packet(TDWP)[8]and quasi-classical trajectory(QCT),[9]was generally based on the accurate potential energy surface(PES).Takayanagi et al.[10]calculated the PESs for the N+H2reaction by using the ab initio multireference configuration interaction(MRCI)method.Pascual et al.[11]employed the modified GROW program developed in Refs.[12]-[15]to select the points and calculated the PES at the MCQDPT2//FORS-MCSCF(7,6)/6-311++G**level of theory by using GAMESS.Later,Adam et al.[6]and Qu et al.[7]proposed the global PESs for4A′′and2A′′states,which are based on the internally contracted MRCI method.Recently,Zhai and Han[16]constructed a high-quality global PES(hereafter denoted as ZH PES)for the ground4A′′state by using the many-body expansion method and the neural network method.For the ZH PES used in this work,all ab initio calculations are carried out in the MOLPRO package and 8645 ab initio energy points are used for fitting the PES with a basis set of aug-cc-pV5Z.The reaction heat of the ZH PES is closer to the experimental value and the predicted potential barrier height is lower than the theoretical value reported in previous work,which can therefore increase the reaction probability,especially at the lower collision energy.

With this ZH PES,the relevant reactions including their isotopic reactions have been studied.Yu X and Yu Y J[17]and Zhang et al.[18]used the QCT method to calculate the polarization-dependent differential cross sections(DCSs)and the three angular distributions of P(θr),P(φr),P(θr,φr)for the reaction N+H2→NH+H.The state-to-state dynamics of the N(4S)+H2(X1Σ+)reaction in ZH PES has been reported by Zhang et al.[19]Furthermore,the reaction H+NH and its isotopic variants have been studied in previous work.Yao et al.[20]studied the state-to-state quantum dynamics of H+NH reaction,calculated and discussed the in fluences of collision energy on the product state-resolved integral cross sections(ICSs)and DCSs.Besides,Li et al.[21]implemented the TDWP and QCT methods to investigate the H+ND reaction and its isotopic reaction.

However,for the title reaction,there are as yet no relevant calculations at the state-to-state level.In this work,we discuss the state-to-state dynamics calculations of the D+ND→N+D2reaction on the4A′′PES.To obtain the more precise information,we calculate state-resolved ICSs and DCSs of the reaction with considering all J values.The rest of the present article is arranged as follows.In Section 2 we brie fly describe the theoretical methods used in our quantum dynamical calculation.Our results are presented and discussed in Section 3.Finally,the main conclusions are dawn from the present study in Section 4.

2.Theoretical method

2.1.PES

In this work,we use an adiabatic ZH PES constructed by Zhai and Han[16]and the PES is a fully dimensional analytic PES composed of 8645 ab initio energy points.The main features of this ZH PES are shown in Fig.1.It can be seen from this figure that there is an energy barrier located at the entrance channel,no well exists and the minimum reaction path occurs in the collinear configuration of the three nuclei.Besides,the title reaction is an exothermic reaction with a heat release of 1.17 eV.More information about ZH PES can be found in Ref.[16].

Fig.1.Potential energy versus reaction path for H+NH reaction on ZH PES.

2.2.Method

When a collision of two molecules happens,its state-tostate quantum dynamics can provide more detailed information and profound observation of the chemical reaction process.But for general systems,the state-to-state quantum scattering calculation takes a long time to obtain the information about accurate dynamics.[22-28]To solve this problem,Zhang and Han[29-31]carefully studied the algorithm of the TDWP method and found that it had a high degree of parallelism.Based on this feature of the algorithm,they modified the program of TDWP into a parallel accelerated version of graphic processing units(GPUs). The relevant calculation shows that this version of GPUs can give high-efficiency dynamical results.[29]In this work,we used the GPUs to perform the calculation about the dynamics of reaction D+ND→N+D2.For the given J value,the Hamiltonian in body- fixed(BF)product Jacobi coordinate can be written as

where R is the distance between N and the center of mass of D2;r is the internuclear distance of diatomic molecule D2;two reduced massμRand μrcan be expressed asμR=2mNmD/(mN+2mD)andμr=mD/2,respectively;is the total angular moment operator;is the rotational angular moment operator of D2;Vpesis the PES used in this work.

The reactant was prepared in the space- fixed(SF)Jacobi coordinate system by using the Gaussian wave packet in R direction and can be expressed as

with

where L=J-j,V=Vpes-Vr(r),Vr(r)is the potential of the product molecule D2.

The radial component of the product wave packet is chosen to be a delta function multiplied by the outgoing asymptotic radial function

where R∞is a fixed radial coordinate of asymptotic region,ν′and j′are the final vibrational and rotational quantum number,respectively.Using this method,the scattering matrix can be obtained in BF Jacobi coordinate as follows:

The coefficients are given by

where i and f denote the initial and final states,h(1)and h(2)are the first and second kind of spherical Hankel functions,respectively.Finally,the state-to-state DCS and ICS can be calculated by the scattering matrix summed over all useful total angular momentum J values,and expressed as follows:

Here,θ is the scattering angle between the reactants and the products,(θ)the reduced rotation matrix,K0and K′are the initial and final projection of the total angular momentum J,andis the transition matrix from initial state ν0j0to final state ν′j′.

3.Result and discussion

3.1.Number aspects

In this work,we perform the state-to-state dynamics calculation for the reaction D+ND→N+D2with the TDWP method developed on the GPUs.To obtain the converged ICSs,the maximum value of J is 49 for collision energies up to 0.5 eV.In order to obtain the result of convergence,many tests have been carried out for ν=0 and j=0 to determine the optimal parameters.The test results are shown in Table 1.It is noted that K is one of the main convergence parameters,which sets an upper limit on the helicity quantum number.From Table 2,we can know that for the selected J and collision energy Ec=0.5 eV,the number of K(K=min(14,J+1))is enough to obtain the converged results.

Table 1.Numerical parameters used in present quantum wavepacket calculations(all parameters are given in a.u.(atomic unit)unless otherwise stated).

Table 2.Convergence test results of title reaction probability as a function of K(with E c=0.5 eV and J=20,30,and 40).

3.2.Total reaction probability and ICS

First,the reaction probability of the title reaction at five selected total angular momenta(J=0,10,20,30,40)are shown in Fig.2.For J=0,because of the existence of the tunneling effect,the threshold energy of the reaction is about 0.03 eV,lower than the barrier height(Ec≈0.08 eV).A broad resonance structure is displayed and it gradually disappears with the increase of J,which has been also found in H+H2[32]and H+Li2[33]reactions.This phenomenon can be partly attributed to the shift of the minimum of PES to higher energy as J increases.Moreover,because of the existence of centrifugal potential,the threshold energy also increases as the J value increases.The total ICSs are displayed in Fig.3.Obviously,the curves of the total ICS are consistent with the results in Ref.[21],so we can conclude that our calculations are reasonable.

Fig.2.Total reaction probabilities versus collision energy for the D+ND→N+D2 reaction at several values of J.

Fig.3.Total ICSs versus collision energy for different product vibrational state distributions in D+ND(ν=0,j=0)reaction.

Fig.4. Plots of(2J+1)weighted opacity function versus J for D+ND→N+D2 reaction at different collision energies.

Figure 4 shows the plots of opacity function versus J for the reaction D+ND at five different values of collision energy.The contribution of J-dependent partial wave to the ICSs is through a 2J+1 factor,which is analogous to the opacity function P(b).As we can see from Fig.4,all the curves show the arch shapes approximately,i.e.,each of the distribution curves reaches the maximum at a certain J value and then decreases to zero at the last accessible value of J.As the collision energy increases,the maximum value of J dependence shifts to the larger values of J.The reactions occur at the smaller values of J at the lower collision energy,and the larger values of J take part in the reaction with the increasing collision energy.Different J dependence corresponds to the different impact parameter b dependence.So,we can think that there are two mechanisms for the title reaction,i.e.,the rebound collision and the stripping collision.The rebound collision is associated with small impact parameter b and leads to the backward scattering,while the stripping collision corresponds to the big impact parameter b and tends to generate the sideward and forward scattering.Overall,we can think that the backward scattering is dominated at low collision energy,and the sideward and forward scattering begin to appear as the collision energy increases.

3.3.Product state distribution

The D2vibrational state-resolved ICSs are also shown in Fig.3.As we can see,ν′=0,1,2,3 vibrationallevels are open in the energy range considered and the D2product is mainly formed in the ν′=0 and ν′=1 states.The curve of ν′=0 increases monotonically as the collision energy increases,while the ν′=1 increases and then decreases with the increase of collision energy.Besides,the result of ν′=1 state is higher than that of ν′=0 before Ec≈ 0.3 eV and the opposite behavior appears after about 0.3 eV.This phenomenon can be due to the fact that the direct mechanism proceeds through a rebound mechanism in which the head-on collision plays an important role at low collision energy.In the head-on collision,the vibrational excitation of the reactant molecule is produced and then transferred to the product molecule.But,the stripping mechanism emerges and dominates as the collision energy increases,and the rotational excitation is more advantageous in this case.

Fig.5.Product rotational state distributions for D+ND(ν=0,j=0)reaction at three different values of collision energy.

Fig.6.(a)3D product vibrational state-resolved reaction cross section and(b)3D product rotational state-resolved cross section.

Figure 5 shows the product rotational state-resolved ICSs at three values of collision energy.All curves reveal the similar features and the peaks lie between j′=0 and the maximum j′value allowed by energy.Unlike the vibrational state distribution,the rotational inversion population appears in all rotational states for the considered energy values.Although the behaviors are similar under different energy values,there are also some differences among them.With the increase of collision energy,the peak value of ICS tends to the higher j′,and the magnitude of the maximum value increases.This indicates that the number of product rotational channels increases as the collision energy increases,and more rotating channels are effectively opened at the higher collision energy.In addition,there is a negative correlation between vibrational excitation and rotational excitation,i.e.,with the increase of product vibrational number ν′,the density of the rotational state decreases gradually.To give a much more indicative view of the product rovibrational state distributions,we plot threedimensional diagram in Fig.6.As shown in this figure,the products mainly concentrate on ν′=0,1 and 0 ＜ j′＜ 15,which is consistent with the distributions shown in Figs.3 and 5.

3.4.DCS

Figure 7 shows the total DCSs of the title reaction at five values of collision energy.It can be seen that the distribution is dominated by the backward scattering at low collision energy.The sideward scattering and the forward scattering begin to appear as the collision energy increases.Besides,the sideward scattering and the forward scattering are obviously enhanced as the collision energy increases.The details of the distribution from the scattering angle 0°to 12°are shown in Fig.7,in which we can find that there is a weak forward scattering at Ec=0.5 eV.This is due to the collision with larger impact parameters involved in the reaction,which fits our guess in Fig.4.From this drawing of partial enlargement,we can find that the DCSs decrease with the decrease of the collision energy in the forward scattering region.When the collision energy is 0.2 eV,the DCSs in the angle range from 0°to 12°are close to zero.As the collision energy continues to decrease to a smaller energy(such as 0.1 eV),the DCS oscillates slightly in the forward scattering region.This can be attributed to the fact that the scattering direction of the product is weakly affected by the smaller collision energy,and the scattering in each direction will have a certain probability.Due to these reasons,the DCS of Ec=0.1 eV is larger than that of Ec=0.2 eV in the small scattering angle region.

Fig.7.Total DCSs versus scattering angle for D+ND→N+D2 reaction at five different values of collision energy.

To show more details of the reaction,we calculate the state-resolved DCSs at three values of collision energy,which are shown in Fig.8.Forlow collision energy case,Ec=0.1 eV,the product experiences predominately the backward scattering in two vibrational states.For the higher energy case,the distribution of the DCS at Ec=0.3 eV is similar to that at Ec=0.5 eV.We can find that the backward scattering relates to the relatively small rotational states.Besides,the sideward scattering is enhanced with the increase of rotational quantum number.

Fig.8.State-resolved DCSs for title reaction at three typical collision energy values.

4.Conclusions

In this work,the efficient GPU version of the TDWP code is used to study the reaction D+ND on a new ab initio PES.The total ICS increases as the collision energy increases after the threshold value,barely having the oscillatory structure.In addition,there is no obvious inversion of quantum number in the vibrational state-resolved ICS,while the rotational state-resolved ICS has inversion phenomenon.Because of the low impact parameter b collisions at low collision energy,the products experience mainly backward scattering.With the increase of collision energy,the sideward scattering and the forward scattering begin to appear and increase gradually,and the product molecules are excited to high rovibrational states.The backward scattering of the product is also observed at high energy,which means that there exist rebound and stripping mechanisms at high energy.