New diabatic potential energy surfaces of the NaH2 system and dynamics studies for the Na(3p) + H2 → NaH + H reaction

The Na(3p) + H2(X1Σg+) → NaH(X1Σ+) + H(2S) reaction plays an important role in the field of diabatic reaction dynamics. A set of new diabatic potential energy surfaces (PESs) of the NaH2 system are structured, which include the diabatic coupling between the lowest two adiabatic states. The electronic structure calculations are performed on the multi-reference configuration interaction level with the cc-pwCVQZ and aug-cc-PVQZ basis sets for Na and H atoms. 32402 geometries are chosen to construct the diabatic data by a unitary transformation based on the molecular property method. The diabatic matrix elements of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{V}}}_{{\bf{11}}}^{{\boldsymbol{d}}}$$\end{document}V11d, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{V}}}_{{\bf{22}}}^{{\boldsymbol{d}}}$$\end{document}V22d and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{V}}}_{{\bf{12}}}^{{\boldsymbol{d}}}$$\end{document}V12d (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{V}}}_{{\bf{21}}}^{{\boldsymbol{d}}}$$\end{document}V21d) are fitted by the artificial neural network model. The spectroscopic constants of diatoms obtained from the present PESs are consistent with the experimental data. The topographical and intersection characteristics of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{V}}}_{{\bf{11}}}^{{\boldsymbol{d}}}$$\end{document}V11d and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{V}}}_{{\bf{22}}}^{{\boldsymbol{d}}}$$\end{document}V22d surfaces are discussed. Based on the new PESs, the time-dependent quantum wave packet calculations are carried out to study the reaction mechanism of the title reaction in detail.

The interactions between electronically excited alkali atoms and hydrogen molecule, including both no-reactive quenching and chemical reactions, have become an interesting topic due to the special advantage for understanding diabatic processes. For the collisions of excited alkali atoms with H 2 , the diabatic couplings necessarily occur, and the conical intersections arise to connect two or more surfaces. Among the issues of reactive collision for these alkali elements, the reactions of K 1-4 , Rb 5,6 , and Cs 7-9 with H 2 proceed by a harpoon model, whereas the reactions of Li [10][11][12][13] and Na [14][15][16] with H 2 follow an insertion mechanism.
The NaH 2 presents a simple prototype of such collision systems, and numerous experimental results [15][16][17][18][19][20][21][22] about the collisions between excited sodium atom and hydrogen molecule are available. Botschwina et al. 17 performed crossed molecular beam experiments on the Na(3p) + H 2 → Na(3 s) + H 2 process. The results of quenching cross section and vibrational energy distribution were presented in their work. Bilililign et al. 15 studied the Na(4p) + H 2 → NaH + H photochemistry reaction by using the "half-collision" pump-probe technique. They observed a bimodal rotational distribution of the NaH product: the minor peak at lower rotational states is due to the repulsive interaction, while the major peak at higher rotational states is attributed to the attractive interaction. Motzkus et al. 20 applied three different nonlinear optical techniques, including coherent anti-Stokes Raman scattering (CARS), resonance-enhanced CARS, and degenerate four-wave mixing, to compare the Na(4p) + H 2 and Na(3p) + H 2 reaction processes. This experiment showed the formation of NaH via the Na(3p) + H 2 reaction follows a two-step mechanism, which is opposite to the direct reaction of Na(4p) + H 2. In 2008, Chang and co-workers 22 studied the rotational and vibrational state distributions of NaH in the reactions of high exited states Na(4 2 S, 3 2 D and 6 2 S) with H 2 by using the pump-probe technique. The authors concluded that the Na(6 2 S) reaction has a dramatically reduced ionization energy, and the corresponding reaction pathway maybe prefer a harpoon model via a near collinear configuration. For the Na(6 2 S) + H 2 reaction, the valence electron of Na hopping mechanism is involved to form an ion-pair Na + H 2 − intermediate. Extensive theoretical studies [23][24][25][26][27] have also been performed on the NaH 2 system based on several potential energy surfaces (PESs) 17,[28][29][30][31][32][33] , which were concentrated on studying the effect of conical intersection and the processes of electronic-to-rovibrational energy transfer. In 1982, Donald et al. 28 calculated the lowest three PESs and diabatic coupling for the Na(3p) + H 2 quenching process by using the diatomics-in-molecules (DIM) method. Then they used a new parametrization for the diabatic potential energy curves of NaH 29 to optimize the preceding PESs, which are only suitable for studying the non-reactive quenching of the Na(3p) + H 2. In 1999, Michael et al. 30 structured an analytical potential energy matrix of the NaH 2 system for the lowest two states based on MP2 calculations, which only include a small number of ab initio energy points near the region of conical intersection. The potential energy matrix can be applied in reaction dynamics calculations of the diabatic processes.
The reaction of lowest electronically excited state Na(3p) with H 2 plays an important role for studying diabatic reaction dynamics correlated with two adiabatic states. The reaction starts at the adiabatic 2 2 A′ surface and then intersects the adiabatic 1 2 A′ surface to enter into the product channel. The formation of NaH molecule by this reaction follows a two-steps process. The first step is the quenching between Na(3p) and H 2 , 2 2 and then Na(3p) collides with the vibrationally excited H 2 , The correctness of this process has been proved by rate equation model 21 . However, the quantum dynamics study of the Na(3p) + H 2 (X 1 Σ g + ) → NaH(X 1 Σ + ) + H( 2 S) reaction has not been reported. For the title reaction, the early diabatic PESs may be not accurate enough to investigate the state-to-state reaction dynamics due to the limit of computational conditions. Fortunately, the recent advances in ab initio theory and the neural network (NN) model make it possible to obtain accurate global PESs. In this work, a set of new diabatic PESs involved two lowest adiabatic states (1 2 A′ and 2 2 A′) and the coupling between them are structured with the NN method. To guarantee the accuracy of the new diabatic PESs, numerous high precision single point energies in a wide coordinate range are calculated, which are used to construct the energy matrix in the diabatic representation by the molecular property method. Then the time-dependent wave packet (TDWP) calculations are carried out on the new diabatic PESs to obtain the quantum dynamics information of the title reaction.

Results
Topographical features of the PESs. Figure 1 presents potential energy curves in the adiabatic and diabatic representations as a function of R Na-HH at four fixed internuclear distances of HH (r HH = 2.1, 2.5, 3.0, 3.5 bohr) in C 2v geometry. The electronic symmetry of the ground and first excited states turn into 1 2 A 1 and 1 2 B 2 in C 2v geometry, respectively. It can be seen that the adiabatic potentials strongly avoid, while the curves in the diabatic representation cross over each other smoothly near the crossing point. The crossing point located at larger R Na-HH with the r HH increasing. Moreover, the adiabatic and daibatic potential energy curves are overlapping when R Na-HH is far away from the crossing point, implying the electronic coupling is zero in the asymptotic region.   Fig. 4. The three-dimensional plots show the fitted PESs are smooth over the whole coordinate space. There is a valley corresponding to Na (3 s) + H 2 channel on the V d 11. For the V d 22 surface, two valleys can be found on the left and right, which represent the Na(3p) + H 2 and NaH + H channels, respectively. Thus in the diabatic representation, the V d 22 surface provides a direct path from Na(3p) + H 2 to NaH + H. For the lower panel of Fig. 4, the red line depicts the position of intersection between the V d 11 and V d 22 PESs. The difference between two diabatic PESs increases in the region away from the intersection position, and this feature is consistent with the results of Fig. 1. The minimum reaction paths for the title reaction at three Na-H-H angles (60°, 90° and 180°) on the V d 22 PES are presented in Fig. 5. For the Na-H-H angle of 60° and 90°, there exist a well along the reaction path, which corresponds the complex-forming reaction. For the Na-H-H angle of 180°, there exists a 0.1 eV height barrier, which implies it is very difficult to initiate the title reaction through the Na-H-H collinear path. Namely, the reaction of Na(3p) with H 2 is dominated by an insertion approach. Combined with previous studies, it could conclude that the reactions of low-lying electronic states Na with H 2 follow an insertion model, whereas the reactions of highly-excited states Na with H 2 could favor the harpoon-type mechanism. Moreover, taking into account the zero point energies of H 2 (0.271 eV) and NaH (0.072 eV) molecules, the endothermicity for forming the NaH molecule is about 0.63 eV.   in the product channel. The threshold increases at large J value due to the increasing centrifugal barrier. Some sharp resonance peaks can be found due to the potential well on V d 22 PES, which gives rise to the formation of intermediate complex. Moreover, the reaction probabilities decrease and the resonance structures become less pronounced as J value increasing, which is attributed to the large centrifugation potential leads to more Na atoms entrance the product channel without the well, thus the lifetime of complex becomes shorter.
In the TDWP calculations, the maximum J value is calculated up to 70, which can ensure the convergence of integral cross sections (ICSs) and differential cross sections (DCSs) when the collision energy is below 1.5 eV. The ICS results of the Na(3p) + H 2 reaction are similar to the Li + H 2 reaction 34 . The ICSs curves are relatively smooth, compared with the oscillating reaction probability curves. The total ICS and several vibrational states (v′ = 0 − 6) ICSs of the title reaction are displayed in Fig. 7. The total ICS value steeply rises at the selected energy region, and the six vibrational excitation channels of the NaH molecule are opened successively. The ICSs of vibrational excitation states keep growth at the collision energy below 1.5 eV. The ground vibrational state ICS rises up to the collision energy reaches 0.94 eV, and then decreases gradually at the energy region from 0.94 to 1.5 eV. It implies more energy are transformed into the internal energy of the NaH molecule, and the product are excited to higher vibrational states.
Furthermore, the product rotational distribution for three selected collision energies (0.8, 1.0, 1.4 eV) are calculated, which are shown in Fig. 8. It can be seen that the rotational states of the product NaH molecule are inverted in all vibrational levels, and the range of rotational quantum number of the product NaH molecule becomes broader with increase of the collision energy. Since more energies can be transferred into the internal energy of the NaH molecule. In all cases, as vibrational quantum increases, the peak of rotational states distribution shifts to lower rotational quantum number. This is because the total energy is constant, and the internal energy shifts from rotation to vibration with increasing vibrational level.
The DCS gives the product angular distribution of a reaction. Figure 9 shows the DCSs of the title reaction at three collision energies (0.8, 1.0, 1.4 eV). It is clear that the product NaH molecule tends to be forward scattered, which implies that the product NaH molecule prefers moving toward the initial direction of the reactant Na atom. This forward bias means that the reaction is dominated by the formation of short-lived complex. With the collision energy increasing, the forward scattered becomes more obvious due to the proportion of direct reactive mechanism increases at a high collision energy.

Discussion
In the present work, a set of new global PESs for the NaH 2 system are constructed in the diabatic representation, which are correlated with the adiabatic 1 A′state and 2 A′state. The ab initio calculations are performed at the internally contracted multi-reference configuration interaction level with the Davidson correction (icMRCI + Q). The aug-cc-PVQZ and cc-pWCVQZ and basis sets are adopted for H atom and Na atom, respectively. The diabatic matrix elements are generated by the transformation of ab initio data based on the molecular property method. A mass of geometries (32402) in a large coordinate space are selected to fit the diabatic PESs using the NN model, and the root mean squared errors (RMSEs) for the diabatic terms  ) → NaH(X 1 Σ + ) + H( 2 S) reaction are carried out to obtain the rigorous quantum dynamics information. The dynamics results show the reaction threshold is consistent with the endothermicity obtained from the diabatic PESs. There exist some oscillation structures on the reaction probability curves due to the complex forming in the potential well. The total ICS steeply rises when the collision energy below 1.5 eV, and the rovibrational states ICSs of the product NaH molecule are presened. In addition, the product NaH molecule tends to be forward scattered, and the forward bias becomes more obvious with increase of the collision energy. As we know, there is no available experimental study which can directly examine the present results. We anticipate our studies could serve as a reference of future experiments for the title reaction.     [35][36][37][38][39][40][41][42] have been developed, and the most important step is to obtain the transformation matrix between the diabatic and adiabatic representations. An effective transformation approach with less calculation burden is to use appropriate molecular properties associated with the transition of electronic states to construct the matrix, and the dipole moment is selected for the NaH 2 system. A brief description about the diabatization scheme is presented below. In this work, the diabatic coupling includes two electronic states, thus the unitary transformation from the adiabatic basis ψ i a to the diabatic basis φ i d can be expressed as where α is called the mixing angle. The energy matrix elements at the diabatic representation is calculated as follows       a  a  a  d  a  d  3  2  3  1  3  2 where ψ a 3 is the adiabatic state 1 2 A″ which is not involved in the electronic coupling, and P denotes the dipole moment operator. The assumption is to make not just for the collinearity. Thus, the α is calculated byα a a a a 3 1

2
Ab initio calculations. The electronic structures of the 1 2 A′ and 2 2 A′ adiabatic states for the NaH 2 system are calculated in C s symmetry using the icMRCI + Q method with a complete active space self-consistent field (CASSCF) reference wave function. In the state-averaged CASSCF calculations, three electronic states (1 A′, 2 A′ and 1 A″) of NaH 2 are assigned to equal weight. Three valence electrons are in fifteen active orbitals (11a′ + 4a′′), consisting of two 1 s orbitals on H, and 3 s, 3p, 3d, 4 s and 4p orbitals on Na. The frozen-core approach is used in the MRCI calculations, and the orbitals (4a′ + 1a′′) are set to close, thus three electrons are included in the calculations of correlation energy. The aug-cc-PVQZ and cc-pWCVQZ basis sets are used for H atom and Na atom, respectively. A total of 32402 geometries are selected to structure the diabatic energy matrix in the Jacobi coordinates. The energy points are defined by 0.6 ≤ r HH /a 0 ≤ 28, 0 ≤ R Na-HH /a 0 ≤ 35, 0 ≤ θ/degree ≤ 90 for the Na−HH region, and 0.8 ≤ r NaH /a 0 ≤ 28, 0 ≤ R H-NaH /a 0 ≤ 35, 0 ≤ θ/degree ≤ 180 for the H−NaH region. All of the ab initio calculations are implemented by MOLPRO package 43 .
Fitting of the diabatic PESs. The energies in the diabatic representation are constructed by combining the adiabatic data with the diabatic transformation based on molecular property method. All of the diabatic energies are used to fit the diabatic terms of V d 11 , V d 22 and V V ( ) d d 12 21 by the NN method. The NN method is an excellent tool to accurately establish PES, and has been used to numerous reactive systems [44][45][46][47][48][49][50] . The feed-forward NN is employed in this work. The NN consists of a set of input signal {x i }, one output signal corresponded to energy and two hidden layers. The output signal of a neuron can be presented as where w i and b j are the weights and bias of interconnecting neurons, respectively. The linear function as the transfer function f(x) in the output layer, and the hyperbolic tangent function is chosen in the two hidden layers, which is written as x x x x The permutation invariant polynomials 51,52 are used in the fitting of each diabatic term. To prevent overfitting, all of the diabatic energies are randomly divided into the training (90%), validation (5%) and testing (5%) sets. The final fitting energy can be expressed as where the superscripts of (1) and (2)  Quantum Dynamics. The quantum dynamics simulation for the Na(3p) + H 2 (X 1 Σ g + ) → NaH(X 1 Σ + ) + H( 2 S) reaction is carried out by the TDWP method on the new NN diabatic PESs, which has been described in detail previously 56,57 . This method is effective to treat the diabatic transition between two electronic states, and only an online involved the main equations is presented below. In the body fixed representation, the Hamiltonian of the Na + H 2 reaction is written as follows  where R is the distance between Na and HH center of mass, and r represents HH bond length. µ R and µ r denote the reduced masses relevant to R and r. J and j are the total and H 2 molecular angular momentums. V is the diabatic PESs of the NaH 2 system. The reactant coordinate based method is used to obtain the S-matrix at the product channel, which is developed by Sun et al. 58 . The states resolved reaction probability can be calculated by ∑ = + .  In this work, only the ground rovibrational state reaction of Na(3p) + H 2 (v 0 = 0, j 0 = 0) → NaH + H are performed the TDWP calculations, and the main dynamics parameters are given in Table 2