In recent decades, the interactions between beryllium atom and hydrogen molecules have been of great attention because of their significance in astrophysics, hydrogen storage, quantum chemistry and other fields. On the one hand, the collision product BeH₂ molecule presents the fundamental and technological interest in potential applications, such as the nuclear materials and rocket fuel technology [1, 2], owing to its small mass and large hydrogen-to-metal mass ratios. Moreover, the molecular BeH₂, with a simple electronic structure, has become an excellent candidate for testing new computational methods for quantum chemistry [3–6]. On the other hand, the further product BeH molecule in the collision process of Be + H₂ is a popular testing target for the electronic structure calculations in open-shell systems [7, 8]. In addition, BeH is also an important interstellar molecule, which has been identified in stars and comets [9, 10].

Various experimental studies on the BeH₂ system have been implemented [11–18]. Tague and Andrews first detected the BeH₂ in molecular form by using infrared spectroscopy and the matrix isolation technique [11]. In their experiment, the pulsed laser evaporated Be atoms react with the hydrogen, and the primary product BeH and BeH₂ are largely favored compared with the other four more complex product molecules of Be₂H, HBeHBeH, HBe(H)₂BeH and HBeBeH. [12] synthesized the gaseous BeH₂ molecule using an electrical discharge facility, which is verified by infrared emission spectroscopy. Their study concluded that the stable BeH₂ is a linearly symmetric molecule with the BeH bond length of 1.334 Å. The high-precision infrared emission spectra of the BeH₂ and BeD₂ molecules were measured by [14]. The antisymmetric stretching modes and some hot bands of the two molecules were studied and the spectroscopic data were accurately determined. In their later study [17], the new vibration-rotation hot bands of the BeH₂ molecule were analyzed, and an accurate value was obtained for the frequency of the bending vibrational mode.

In the theoretical aspect, numerous ab initio calculations on the BeH₂ molecule have also been performed [19–26]. Martin and Lee [19] accurately calculated the quartic force field of BeH₂ using the CCSD(T) method, and the obtained spectroscopic constants are consistent with the corresponding experimental measurements. Hrenar et al. reported the first potential energy surface (PES) of the BeH₂ molecule used for the vibrational configuration-interaction calculations by a multilevel scheme [23]. Their calculated results can reproduce the experimental values of the gas phase measurements and matrix isolation. The ground-state equilibrium structure and PES of BeH₂ were calculated utilizing the CCSD(T) method combined with the cc-pVTZ through cc-pV6Z basis sets by Koput and Peterson [24]. Furthermore, the rovibrational energy levels of BeH₂ and its isotopic variations of BeD₂ and BeHD were accurately calculated by a variational method. The newest PES of the BeH₂ system was constructed by Li and Roy [25] utilizing the three-dimensional spline interpolation over 6,864 energy points with the internally contracted multi-reference configuration interaction (icMRCI)/aug-cc-pV5Z level. On this PES, the spectral constants of the BeH₂ and BeD₂ molecules were accurately calculated and the corresponding data of the BeHD molecule were predicted.

Although the BeH₂ system has received great attention both experimentally and theoretically, most of those studies focused on its structural and spectral properties, and the dynamics mechanisms of the Be + H₂ reaction process have not been reported up to now. In theory, the most reliable approach for obtaining the accurate dynamics information of a chemical reaction is to implement rigorous quantum scattering calculations on a globally high-precision PES [27, 28]. The previous PESs of the BeH₂ system are extremely reliable and accurate for describing the BeH₂ complex, whereas they are not suitable for the reaction dynamics calculations since some key regions where the reaction could reach are not included. Therefore, constructing a global and accurate BeH₂ PES is a crucial premise for studying the microscopic dynamics mechanisms of this reactive system.

Herein, a high-fidelity ground-state BeH₂ PES is represented based on a mass of high-precision ab initio energy points and the permutation invariant polynomial-neural network (PIP-NN) scheme [29, 30]. Moreover, the quantum dynamics calculations at the state-resolved level for the Be(¹S) + H₂(v ₀ = 0, j ₀ = 0) → BeH + H reaction are carried out by the time-dependent wave packet (TDWP) method [31, 32] on this newly constructed PES. The computational details and the characteristics of the PES are given in Section 2. Section 3 displays the calculated dynamics results and the relevant discussion of the dynamics mechanisms for the title reaction and Section 4 concludes this work.

For most of the triatomic and some tetratomic reactive systems, the quantum TDWP method [31, 32, 53–55] is a high-efficiency and accurate tool for calculating the dynamics data. The full-dimensional quantum dynamics calculations of the Be(¹S) + H₂(v ₀ = 0, j ₀ = 0) → BeH + H reaction are carried out on this newly constructed PIP-NN PES by the TDWP method for understanding the state-resolved dynamics mechanisms. The Coriolis coupling effect is included in the quantum TDWP calculations. Here, only the main equations in the TDWP calculations are displayed below. The Hamiltonian of the title reaction can be expressed as: H ^ = − ℏ 2 2 μ R ∂ 2 ∂ R 2 − ℏ 2 2 μ r ∂ 2 ∂ r 2 + ( J ^ − j ^ ) 2 2 μ R R 2 + j ^ 2 2 μ r r 2 + V ^ (9) where µ _r and µ _R are the reduced masses associated with r and R in the Jacobi coordinate, respectively. J and j express the total angular momentum quantum number of BeH₂ and rotational angular momentum quantum number of H₂, respectively. The initial wave packet consists of a Gaussian type wave function, a rovibrational eigenfunction of H₂, and an eigenfunction of the total angular momentum, written as: Ψ a v 0 j 0 l 0 J M ε ( t = 0 ) = G ( R α ) ϕ v 0 ( r α ) j 0 ( θ α ) | J M j 0 l 0 ε 〉 (10)

To avoid the reflection of wave packet at the grid edge, the absorption potential used in the TDWP calculations is defined as: D ( x ) = { exp [ − Δ t ⋅ C a ⋅ ( x − x a x b − x a ) 2 ] , x a ≤ x ≤ x b exp [ − Δ t ⋅ C b ⋅ ( x − x b x e n d − x b ) 2 ] × exp ( − Δ t ⋅ C a ) , x b < x ≤ x e n d (11) where C _i and x _i (i = a, b) represent the strength and positions of the absorption potential, respectively. Here, the time evolution of the wave packet is realized by the split operator scheme [56] and using the reactant coordinate-based method [57, 58] to extract the state-resolved S-matrix. The rovibrationally state-resolved reaction probability obtained by the S-matrix is expressed as: P v j ← v 0 j 0 J = 1 2 j 0 + 1 ∑ K ∑ K 0 | S v j K ← v 0 j 0 K 0 J | 2 (12)

The state-resolved integral cross sections (ICSs) are calculated by summing the probabilities of all the calculated partial wave J: σ v j ← v 0 j 0 = π ( 2 j 0 + 1 ) k v 0 j 0 2 ∑ K ∑ K 0 ∑ J ( 2 J + 1 ) | S v j K ← v 0 j 0 K 0 J | 2 (13) where k v 0 j 0 is the momenta in the entrance channel. The state-resolved differential cross sections (DCSs) can be obtained by the following equation: d σ v j ← v 0 j 0 ( ϑ , E ) d Ω = 1 ( 2 j 0 + 1 ) ∑ K ∑ K 0 | 1 2 i k v 0 j 0 ∑ J ( 2 J + 1 ) d K K 0 J ( ϑ ) S v j K ← v 0 j 0 K 0 J | 2 (14) where d K K 0 J ( ϑ ) represents the reduced Wigner matrix, and ϑ express the scattering angle.

In this work, the reactant H₂ molecule is set at its ground-rovibrational state of v ₀ = 0, j ₀ = 0, and the number of partial waves is calculated up to 65, which can obtain the convergent ICS and DCS up to the collision energy of 5.0 eV. In Table 3, the main parameters determined by many times tests of convergence in the TDWP calculations are listed.

The collision energy dependence of total reaction probabilities for the Be(¹S) + H₂(v ₀ = 0, j ₀ = 0) → BeH + H reaction with four partial waves (J = 0, 20, 40 and 50) are presented in Figure 7A. For J = 0, the curve exhibits relatively dense oscillation structures, which are attributed to the potential well on the reactive path. The title reaction is dominated by the global MEP and there is a well with the depth of 1.632 eV, resulting in obvious quantum resonances because numerous bound and quasi-bound states can be formed in the well. As the increase of J values, the reactive threshold becomes larger and the oscillations are gradually weakened. This is because the increasing centrifugal barrier reduces and even smooths the effective potential well, and the other collision channels shown in Figure 5 are opened, causing the amplitudes of oscillations on the reaction probability curves become less pronounced. Figure 7B shows the collision energy dependence of total ICS for the title reaction. The total ICS value increases monotonically with the increase of collision energy, which is consistent with the characteristic of an endothermic reaction. Compared to the reaction probabilities, there is no oscillation structures on the ICS curve due to the superposition of all the calculated partial waves.

To understand the dynamics mechanisms of the Be(¹S) + H₂(v ₀ = 0, j ₀ = 0) → BeH + H reaction at the state-to-state level, the rovibrationally state-resolved ICSs of the product BeH molecule at four collision energies (3.0, 4.0, 4.5, and 5.0 eV) are shown in Figure 8. For the collision energy of 3.0 eV, the BeH molecule only can be excited to the lowest three vibrational states, but the maximum of the rotational quantum number can reach j′ = 21 at v′ = 0, and the peak value of the rovibrationally state-resolved ICS is located at v′ = 0, j′ = 16. The presented vibrationally cold and rotationally hot distribution conforms to the complex-forming mechanism. More rovibrational states become available with the increase of collision energy, and there is a population inversion of the vibrational quantum number. For the collision energy of 5.0 eV, the product BeH molecule can populate at very high rovibrational states (v′ = 10, j′ = 28), suggesting more collision energy is effectively transformed into the internal energy of the product molecule. The contributions of high-order partial waves are larger and more reaction paths are gradually opened as the collision energy increases, thus the lifetime of the forming BeH₂ complex becomes shorter and the title reaction prefers a direct H-abstraction process when the collision energy is large enough.

To study the dynamics process of the Be(¹S) + H₂(v ₀ = 0, j ₀ = 0) → BeH + H reaction more intuitively by giving the angular distribution of the product molecule, Figure 9 presents the total DCSs varying with the scattering angle and collision energy. It is clear that the peak values of the angular distribution are located at the two extreme angles (0° and 180°) and the forward-backward symmetric DCSs are displayed when the collision energy is slightly larger than the reactive threshold, which is due to the forming of a BeH₂ complex supported by the potential well on the global MEP. With the increase of collision energy, the product BeH molecule increasingly prefers the forward scattering, showing an obviously non-statistical behavior. It also can be explained by the increasing contributions of the centrifugal barriers and more open reactive paths without a well at large collision energy. The calculated results of the total DCS further imply that the title reaction follows the complex-forming mechanism near the reactive threshold, whereas s direct H-abstraction process gradually plays a dominant role at high collision energy.

In this paper, a globally accurate ground-state BeH₂ PES is structured using the PIP-NN scheme based on 12371 ab initio points calculated at the icMRCI + Q/aug-cc-pCV5Z level. The PES can accurately reproduce the original ab initio data in each region, and the global fitting RMSE is only 1.972 meV. The molecular constants of H₂(X¹Σ_g ⁺) and BeH(X¹Σ⁺) calculated on the PES are consistent with the corresponding experimental data, and the PES can reproduce the characteristics of stationary points well. The GM and TS of the ground-state BeH₂ correspond to the D _∞h and C _2v symmetries, respectively. The topographic features of the PES are described in detail. On this newly constructed PES, the dynamics calculations are performed on the Be(¹S) + H₂(v ₀ = 0, j ₀ = 0) → BeH + H reaction at the state-to-state level by the quantum TDWP method for understanding the microscopic mechanisms. The endothermicity of the title reaction determined by the PES is 2.716 eV. There exist obvious oscillation structures on the curves of reaction probabilities since the well on the global MEP can support numerous bound and quasi-bound states, and the total ICS increases monotonically with the increase of collision energy. The rovibrationally state-resolved ICSs present vibrationally cold and rotationally hot distribution at relatively low collision energy, and the product BeH molecule can populate at very high rovibrational states. The total DCSs are forward-backward symmetric when the collision energy is slightly larger than the reactive threshold, but only the forward scatting is presented at high collision energy. The dynamics results indicate that the title reaction follows the complex-forming mechanism near the reactive threshold, whereas a direct H-abstraction process gradually plays the dominant role at high collision energy. Further dynamics studies for this reaction system can be carried out on the presented PES, such as the effects of rovibrational excitations and isotopic substitutions of the H₂ molecule, and the dynamics data calculated in this paper would be of importance in the experimental studies on the title reaction.