A comparative study of the low energy HD+o-/p-H2 rotational excitation/de-excitation collisions and elastic scattering

The Diep and Johnson (DJ) H$_2$-H$_2$ potential energy surface (PES) obtained from the first principles [P. Diep, K. Johnson, J. Chem. Phys. 113, 3480 (2000); 114, 222 (2000)], has been adjusted through appropriate rotation of the three-dimensional coordinate system and applied to low-temperature ($T<300$ K) HD+$o$-/$p$-H$_2$ collisions of astrophysical interest. A non-reactive quantum mechanical close-coupling method is used to carry out the computation for the total rotational state-to-state cross sections $\sigma_{j_1j_2\rightarrow j'_1j'_2}(\epsilon)$ and corresponding thermal rate coefficients $k_{j_1j_2\rightarrow j'_1j'_2}(T)$. A rather satisfactory agreement has been obtained between our results computed with the modified DJ PES and with the newer H$_4$ PES [A.I. Boothroyd, P.G. Martin, W.J. Keogh, M.J. Peterson, J. Chem. Phys. 116, 666 (2002)], which is also applied in this work. A comparative study with previous results is presented and discussed. Significant differences have been obtained for few specific rotational transitions in the H$_2$/HD molecules between our results and previous calculations. The low temperature data for $k_{j_1j_2\rightarrow j'_1j'_2}(T)$ calculated in this work can be used in a future application such as a new computation of the HD cooling function of primordial gas, which is important in the astrophysics of the early Universe.


I. INTRODUCTION
The simplest quantum-mechanical scattering problem, that between two hydrogen molecules, has been a formidable challenge to chemists and physicists working in this area. This is because the elastic and inelastic H 2 +H 2 scattering have all the complications of a complex molecular scattering process like rotational and vibrational excitations and de-excitations. Many theoretical and experimental methods have been developed and used for the study of these processes [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] . Additionally, there has been great interest in the study of HD+H 2 scattering. Although, this system bears some similarity with H 2 +H 2 , the identical system symmetry is broken here and hence many aspects of the theoretical formulation can be tested in the system under different symmetry conditions 4, [16][17][18][19][20][21][22] .
The PESs of H 2 +H 2 and HD-H 2 are basically the same six-dimensional functions. The fact follows from a general theoretical point of view and the Born-Oppenheimer model 23 . However, from an intuitive physical point of view, the two collisions should have rather different scattering output. This is because the H 2 and HD molecules have fairly different rotational constants and rotationalvibrational spectrum. They have different internal symmetries and dipole moments. Collisional properties of a) Electronic mail: rasultanov@stcloudstate.edu b) Electronic mail: dcguster@stcloudstate.edu c) Electronic mail: adhikari@ift.unesp.br; http://www.ift.unesp.br/users/adhikari these systems are expected to be highly sensitive to dipole moments. Further, the HD-H 2 PES can be derived from H 2 +H 2 by adjusting the coordinate of the center of mass of the HD molecule. Once the symmetry is broken in H 2 +H 2 by replacing the H by the D atom in one H 2 molecule one can obtain the HD-H 2 PES. The new potential has all parts of the full HD-H 2 interaction including HD's dipole moment. However, in a case when a PES is formulated for fixed interatomic coordinates between the two hydrogen atoms in the H 2 -H 2 system, it would be difficult to extract the PES of HD-H 2 , because the position of the HD center of mass is different from that of H 2 . In this work we apply a rotational method for calculation of the HD-H 2 PES from that of H 2 -H 2 . The method is based on a rotation of the three-dimensional space from the body-fixed H 2 -H 2 coordinate system to appropriately adjust the HD molecule center of mass.
The hydrogen molecule plays an important role in many areas of astrophysics. For example, the interstellar medium (ISM) cooling process is associated with the energy loss after inelastic collisions of the particles. These processes convert the kinetic energy of ISM's particles to their internal energies: for instance, in the case of molecules to their internal energy of rotationalvibrational degrees of freedom. Therefore, in order to accurately model the thermal balance and kinetics of ISM one needs accurate state-to-state cross sections and thermal rate constants. Theoretical state-resolved treatment requires precise PES of H 2 -H 2 , and a reliable dynamical method for H 2 +H 2 /HD collision 5,16,21 . However, on the other hand, experimental measurements of these cross sections is a very difficult technical problem. Unfortu-nately, up to now no reliable experiments are available on these collision processes, which have important astrophysical applications. Moreover, different calculations with various H 2 -H 2 PESs showed rather different results for important H 2 +H 2 thermal rate coefficients.
The possible importance of the HD cooling in ISM was first noted in 1972 by Dalgarno and McCray 24 . Because of the special properties of HD this molecule is even more effective cooler than H 2 . This happens at low temperatures T < 100 K 24 , where the cooling function of HD becomes ∼15 times larger than the H 2 cooling function 24 . As we mentioned in the previous paragraph, to carry out calculation of the cooling function one needs to know precise cross sections and thermal rate coefficients of the rotational-vibrational energy transfer collisions. That is why there is a constant interest to reliable quantum-mechanical computation of different rotationalvibrational atom-molecular energy transfer collision cross sections and corresponding thermal rate coefficients [25][26][27][28][29] .
A realistic full-dimensional ab initio PES for the H 2 -H 2 system was first constructed by Schwenke 30 and that potential was widely used in a variety of methods and computation techniques. Flower 31 have produced refined PES for the H 2 -H 2 system. These PESs have been used in several different calculations [32][33][34][35][36] . However, in our previous works 12, 34 we found that in the case of low energy H 2 +H 2 collision these PESs provide different results for some specific state-resolved cross sections. The difference may be up to an order of magnitude. This fact was also confirmed in work 13 . The BMKP PES, probably, needs future improvements 12,13,34 , because the DJ PES gives better agreement with existing experiments for H 2 +H 2 than the BMKP potential.
Nonetheless, as a first trial Sultanov et al. applied the BMKP PES to the low-energy HD+H 2 collisions after appropriate modification 37,38 . When the results of these calculations 37,38 were compared with prior studies 5,21 a relatively good agreement was found with Schaefer's results 5 . At the same time substantial differences with the newer data obtained by Flower was noted 21 . Therefore, in the light of these circumstances, it would be useful to apply another modern H 2 -H 2 PES 1 to HD-H 2 collisions. The DJ potential was formulated specifically for the symmetric H 2 -H 2 system when distances between hydrogen atoms are fixed at a specific equilibrium value in each H 2 molecule. In this paper we provide the first calculation describing collisions of rotationally excited H 2 and HD molecules using the DJ PES. We also provide a comparison with our previous calculations 37, 38 with the BMKP PES, which is a global six-dimensional surface.
In Sec. II we briefly outline the quantum-mechanical close-coupling approach and our method to convert the symmetric DJ PES to be appropriate for calculations of the HD+H 2 system. We represent results for selected state-to-state cross sections and thermal rate co-efficients for rotational excitation/de-excitation of HD in low-energy collisions with o-/p-H 2 in Sec. III. Finally, in Sec. IV we present a brief summary of our findings and conclusion.

A. Dynamical equations
The Schrödinger equation for the (12)+(34) collision in the center-of-mass frame, where 1, 2, 3 and 4 are atoms and (12) and (34) are diatomic molecules modeled by linear rigid rotors, is 3,39,40 where P R3 is the momentum operator of the kinetic energy of collision, R 3 is the collision coordinate, whereas R 1 and R 2 are relative vectors between atoms in the two diatomic molecules as shown in the Fig. 1, and LR 1(2) are the quantum-mechanical rotation operators of the rigid rotors. Here, is the reduced mass of the two diatomic molecules (12) and (34) and µ 1(2) ≡ [m 1(3) m 2(4) ]/[(m 1(3) + m 2(4) )] are reduced masses of the two molecules. The vectorsR 1(2) are the angles of orientation of rotors (12) and (34), respectively; V ( R 1 , R 2 , R 3 ) is the PES for the four-atomic system (1234), and E is the total energy in the center-of-mass system. The linear rigid-rotor model used in this calculation was already applied in some previous studies 21,39,40 . For the considered range of kinetic energies of astrophysical interest the model is quite justified.
The eigenfunctions of the operators LR 1 (2) in Eq. (1) are simple spherical harmonics Y jimi (r). To solve this equation the following angular-momentum expansion is used: where U JM j1j2j12L (R 3 ) are unknown coordinate functions, J is total angular momentum quantum number of the (1234) system and M is its projection onto the space fixed Z axis, and channel expansion functions are the following: The vector R3 connects the center of masses of the HD and H2 molecules and is directed over the axis OZ, θ1 is the angle between R1 and R3, θ2 is the angle between R2 and R3, ϕ1 is the torsional angle, j1, j2, L are quantum angular momenta over the corresponding Jacobi coordinates R1, R2, and R3.
here j 1 + j 2 = j 12 , j 12 + L = J,m 1 ,m 2 ,m 12 andm are projections of j 1 , j 2 , j 12 and L, respectively. The C's are the appropriate Clebsch-Gordan coefficients. The quantum-mechanical momenta j 1 , j 2 and L over corresponding Jacobi vectors R 1 , R 2 and R 3 are also shown in Fig. 1. Upon substitution of Eq. (2) into Eq. (1), one obtains a set of coupled second order differential equations for the unknown radial functions where α ≡ (j 1 j 2 j 12 L). We apply the hybrid modified logderivative-Airy propagator in the general-purpose scattering program MOLSCAT 41 to solve the coupled radial Eq. (4). The log-derivative matrix is propagated to large intermolecular distances R 3 , since all experimentally observable quantum information about the collision is contained in the asymptotic behavior of functions U JM α (R 3 → ∞). The numerical results are matched to the known asymptotic solution to derive the scattering S-matrix S J αα ′ : where k αα ′ = 2M 12 (E + E α − E α ′ ) is the channel wave number, E α(α ′ ) are rotational channel energies for α and α ′ . The method was used for each partial wave until a converged cross section was obtained. It was verified that the results are converged with respect to the number of partial waves as well as the matching radius, R 3max , for all channels included in our calculation. The cross sections for rotational excitation and relaxation can be obtained directly from the S-matrix. In particular the cross sections for excitation from j 1 j 2 → j ′ 1 j ′ 2 summed over the finalm ′ 1m ′ 2 and averaged over the initialm 1m2 are given by The relative kinetic energy of the two molecules in the center of mass frame is: where B 1 (2) are the rotation constants of rigid rotors (12) and (34), respectively, of total angular momentum j 1 (2) . The relationship between the rotational thermal-rate coefficient k j1j2→j ′ 1 j ′ 2 (T ) at temperature T and the corresponding cross section σ j1j2→j ′ 1 j ′ 2 (ǫ), can be obtained through the following weighted average 42 where k B is Boltzman constant and ǫ s is the minimum relative kinetic energy of the two molecules for the levels j 1 and j 2 to become accessible.

B. The BMKP H2-H2 PES
The BMKP PES 31 , is a global six-dimensional potential energy surface for two hydrogen molecules. It was especially constructed to represent the whole interaction region of the chemical reaction dynamics of the fouratomic system and to provide an accurate van der Waals well. In the six-dimensional conformation space of the four-atomic system the PES forms a complicated threedimensional hyper surface 31 . To compute the distances between the four atoms R 1 , R 2 and R 3 , the BMKP PES uses Cartesian coordinates. Therefore it was necessary to convert spherical coordinates used in the closecoupling method 41 to the corresponding Cartesian coordinates and compute the distances between the four atoms followed by calculation of the PES 37,38 . Without the loss of generality the procedure used a specifically oriented coordinate system OXY Z 37 . First we introduce the Jacobi coordinates { R 1 , R 2 , R 3 } and the radius-vectors of all four atoms in the space-fixed coordinate system OXY Z: { r 1 , r 2 , r 3 , r 4 } (not shown in Fig.  1). As in Ref. 37 we apply the following procedure: the center of mass of the HD molecule is put at the origin of the coordinate system OXY Z, and the R 3 is directed to center of mass of the H 2 molecule along the OZ axis, as shown in Fig. 1. Then R 3 = {R 3 , Θ 3 = 0, Φ 3 = 0}, with Θ 3 and Φ 3 the polar and azimuthal angles, Next, without the loss of generality, we can adopt the OXY Z system in such a way, that the HD interatomic vector R 1 lies on the XOZ plane. Then the angle variables of R 1 and R 2 are: One can see, that the Cartesian coordinates of the atoms of the HD molecule are: Defining ζ = m 4 /(m 3 + m 4 ), we have and the corresponding Cartesian coordinates: C. The Modified Diep and Johnson H2-H2 potential Below we briefly present our method to convert the symmetric DJ H 2 -H 2 PES 1 to be suitable for the nonsymmetric system HD+H 2 . The method is based on a mathematical transformation technique, i.e. a geometrical rotation of the three-dimensional (3D) space and the corresponding space-fixed coordinate system OXY Z 43 , as shown in Fig. 2. It is more convenient here to use a slightly different from Fig. 1 orientation of the coordinate system OXY Z: now the center of mass of the H 2 molecule is set at the origin of the space-fixed OXY Z, and R 3 ( R ′ 3 ) is directed to the center of mass of H 2 (HD). The few-body system (1234) can be characterized by four radius-vectors: { r 1 , r 2 , r 3 , r 4 }, or alternatively by so-called three Jacobi vectors: { R 1 , R 2 , R 3 }. Usually the second option is more convenient for describing the quantum-mechanical few-body systems. Next, the initial geometry of the system is taken in such a way that the Jacobi vector R 3 connects the center of masses of the two molecules and is directed over the OZ axis. We can also choose OXY Z in such a manner that the Jacobi vector R 2 lies in the X-Z plane. Finally, the vector R 1 can be directed anywhere. Then the spherical coordinates of the The geometrical configuration of the HD+H2 system in the framework of two different Cartesian coordinate systems: OXY Z and OX ′ Y ′ Z ′ . Here: x is the distance between OH 2 and OHD, η is the rotation angle over the OZ axis, i.e. the angle between the axes OZ and OZ ′ and between the axes OX and OX ′ , r1, r2, r3 and r4 are the radius-vectors of the corresponding atoms. The vectors R2, R3, and R ′ 3 , belong to the XOZ plain.
First we start with the original DJ PES of Eqs. (15)− (19). Then we replace one hydrogen atom "H" with a deuterium atom "D". Thus we break the symmetry by shifting the center of mass of one H 2 molecule to another point, that is from O H2 to O HD as shown in Fig. 2, we also need to take into acount the difference between the number of rotational states in H 2 +H 2 and HD+H 2 . The length of the vector x is x = | R 2 |/6. Now we rotate the OXY Z coordinate system around the OY axis in such a way that the new OZ ′ axis goes through the point O HD . The OY ′ axis and the old OY axis are parallel, and the angle of this small rotation is η, see Fig. 2. This transformation converts the initial Jacobi vectors in OXY Z: As a result of this simple procedure we obtain a new PES, namely: It is quite obvious, that the rotation does not affect the coordinate function V l1,l2,l (R 3 ) in Eq. (15). Note, in a previous calculation of low energy collision between monodeuterated ammonia NH 2 D and a helium atom, a somewhat similar spatial rotational-translational procedure has been applied to the original PES of the NH 3 -He system 44 .
Next, any rotation of the 3D OXY Z coordinate system can be represented by Euler angles, i.e. {α, β, γ} 45,46 . In this work we choose the following Euler angles: {α = 0, β = η, γ = 0}. To calculate the value of the rotational angle η in Fig. 2, one can consider the internal triangle △O HD OO H2 which is shown in Fig. 2. The angle η is determined from the following formula: The derivation of (21) can be expressed as follows. First, the angles of the triangle △O HD O H2 O satisfy the following equation (π − θ 2 ) + η + θ ′ 2 = π or θ 2 = η + θ ′ 2 . Secondly, based on the law of sines for △O HD OO H2 : where ε = π − θ 2 . Because sin ε = cos θ 2 we have: and finally: from which one can directly obtain the expression (21). In such a way the rotation of the coordinate system from OXY Z to O ′ X ′ Y ′ Z ′ makes a corresponding transformation of the coordinates of the atoms in the 4-body system and the distance between two molecules. One has the following relations among old and new variables 45 : In the calculation of HD+H 2 with the DJ PES one has to use new coordinates θ ′ 1 , θ ′ 2 , φ ′ 12 , R ′ 3 . However, the potential (15) has been expressed through the old H 2 -H 2 variables, hence it is to be transformed to new variables using Eqs. (25) − (28).

III. NUMERICAL RESULTS
Results for the low energy HD+o-/p-H 2 elastic scattering cross-sections and few selected quantum-mechanical rotational transitions together with the corresponding results from previous studies 37,38 are presented below. The cross sections are also compared with the corresponding results from 5 . Additionally, we compare our results for the thermal rate coefficients (8) with previous calculations 5, 21 . The results for the low energy elastic scattering cross sections cannot be compared with other theoretical/experimental data. To the best of our knowledge such calculations do not exist.

A. Elastic scattering
In Fig. 3 we show the present elastic scattering cross sections for HD+o-/p-H 2 collisions at low and ultra-low energies calculated using three different potentials: the BMKP potential of Sec. II B, the modified DJ potential of Sec. II C, and the original DJ potential appropriate for the H 2 -H 2 system. In each cross section there is a prominent resonance. In Fig. 3 although the shapes of three cross sections are similar to each other the resonance peak and the position of the resonance differ significantly when we use the original DJ potential without the modifications described in the Sec. II. This is because the original, symmetrical DJ potential 1 does not have all the asymmetrical features of the HD+H 2 interaction. Moreover, these features are of crucial importance for HD+H 2 scattering. In Table III we show the elastic cross sections σ el (ǫ) for a few selected energies calculated with the three potentials. We also include in Table III the results for scattering lengths a 0 's. As can be seen from Fig. 3 and Table III we have obtained extremely good agreement between our calculations with the BMKP and with the modified DJ PESs. For the considered range of the kinetic energies we obtained full numerical convergence, for instance, the total angular momentum J TABLE I. The elastic scattering cross sections σ el (10 −16 cm 2 ) at selected relative kinetic energies ǫ (cm −1 ) and corresponding scattering lengths ascatt (10 −8 cm) in the o-/p-H2 + HD → o-/p-H2 + HD low energy collisions calculated with three different potentials: modified DJ PES from Sec. II C and original BMKP 31  was used up to the maximum value J max = 12 in this calculation.

B. Non-elastic channels
The main goal of the present study is not to obtain results for every possible transition cross section in the HD+H 2 collision, but rather to demonstrate how the appropriately undertaken 3D rotation of the symmetrical surface of the H 2 -H 2 system could be adapted for collision of the H 2 and HD molecules. We carry out calculations for a fairly wide region of the collision energies, i.e. from 3 K to up to 300 K. This temperature interval is relevant for future calculations of the important HDcooling function 24 . Here we compute few transition cross sections in which we noticed substantial differences between our results obtained with the two different PESs from Refs. 1 and 31 and also between our results and the data from previous studies 5, 21 .
In Figs. 4 to 8 we show these results for a few selected state-to-state HD+o-/p-H 2 integral cross sections. The results have been calculated with two different potentials, specifically, with the BMKP PES 31 of Sec. II B and the modified DJ PES of Sec. II C. When possible we compare our results with existing previous calculations for the state-selected total cross sections and thermal rate coefficients. It is necessary to mention that the original DJ surface 1 cannot be correctly used for the non-elastic or transition channel calculations, i.e. for the rotational state transitions in an HD+H 2 collision. If the original, unmodified DJ potential is applied one can get extremely low numbers (∼ 10 −35 ) for the HD+H 2 rotational stateto-state probabilities and cross sections. That is why it was necessary to undertake the geometrical modifications to the DJ surface described in Sec. II C. We obtained a fairly good agreement between our calculations with the use of the two PESs. Besides some qualitative differences in the state-resolved total cross sections the overall behavior of the cross section σ j1j2→j ′ 1 j ′ 2 (v) was found identical, where v is the relative velocity of the two molecules. However we obtained substantial differences between our σ j1j2→j ′ 1 j ′ 2 (v) cross sections and corresponding data from Ref. 5.
In Fig. 4 we plot the total cross sections for the rotational transitions (02) - (20) and (13) - (11). The notation of the rotational quantum numbers can be understood by comparison to the following equations: where the numbers in parenthesis denote rotational quantum numbers j 1 , j 2 etc.. On the upper plot it is seen that while the two results of this study are in fairly good agreement between each other we obtain substantial disagreements with the result of Ref. 5. The same is true in the lower plot, although the behavior of these cross sections has a common character.
In Fig. 5 we show the cross sections for rotational transitions in HD+H 2 for the BMKP and the modified DJ PESs for the processes HD(1) + H 2 (3) → HD(0) + H 2 (1), and HD(1) + H 2 (3) → HD(2) + H 2 (1). The corresponding results from Ref. 5 are also shown. We again obtain a fairly good agreement between our results computed with the BMKP and the modified DJ PESs, however, differing significantly from the corresponding Schaefer result 5 . Unfortunately we cannot compare these cross sections with the results of the calculation by Flower 21 , in which a different H 2 -H 2 potential from 30 was used. This is because Flower's data includes results mostly for the thermal rate coefficients. However, within the next subsection in Table II we compare our few selected rotational state-resolved thermal rate coefficients with the corresponding results from Refs. 5 and 21.
From the astrophysical point of view, perhaps, one needs only precise rotational and in some rare cases vibrational state-to-state thermal rate coefficients k j1j2→j ′ 1 j ′ 2 (T ) in H 2 +H 2 , HD+H 2 etc. These quantities are less sensitive to interaction potentials. However, the overall behavior of all possible state-selected cross sections should be very important for calculation of the thermal rates as seen in Fig. 4. It appears that the cross sections are much more sensitive to the PESs used in the calculations. Hence it is useful and even probably important in some specific cases to compare the cross sections from various calculations where different PESs have been used.

C. Thermal rate coefficients
In Table II our rotational thermal rate coefficients k j1j2→j ′ 1 j2 (T ) are presented. In this paper we choose only three rotational double-transitions in the collision: (02)-(20), (12)- (20), and (20)-(02). These results were obtained from corresponding state-resolved integral cross sections σ j1j2→j ′ 1 j2 (ǫ) with the use of the expression (8). In a previous study 21 , Flower compared his rotational transition thermal rate coefficients with the corresponding Schaefer data 5 and found substantial differences between his results and results from Ref. 5 even at high temperatures. The point is that in these processes the transition occurs in the two molecules simultaneously. In such rotational transitions the probabilities and cross sections should be very sensitive to the interaction potential. Perhaps, because of this reason the results of Refs. 5 and 21 differ so dramatically for transitions like (02) - (20).  In fact, we also obtained substantial differences between our thermal rates and with corresponding results from Refs. 5 and 21. Particularly the large differences were detected at the lower temperature regime.
In Table II we show our thermal rate coefficients k j1j2→j ′ 1 j ′ 2 (T ) for the processes HD(0)+H 2 (2) → HD(2)+ H 2 (0), and HD(1) + H 2 (2) → HD(2) + H 2 (0), and those for HD(2) + H 2 (0) → HD(0) + H 2 (2). We present our results calculated with the BMKP 31 and with the modified DJ 1 PESs. As before, our results computed with these two potentials are close to each other. However, one can see that our results and results from Refs. 5 and 21 differ significantly. This happens particularly at low temperatures, for instance from 10 K to 50 K. Finally in this section, for the process HD(2) + H 2 (0) → HD(0) + H 2 (2) the difference between our k j1j2→j ′  the authors also used two different PESs and different quantum-mechanical methods. However, we would like to point out here that our results computed with two newer PESs 1,31 are relatively close to each other.
D. Application of the detailed balance principle By using the time reversibility (reciprocity) principle one can obtain the detailed balance equation for the direct and reverse energy transfer processes or reactions 47 .
In the case of the inelastic scattering a + b ⇋ a ′ + b ′ the detailed balance principle 47 relates the direct a + b and the reverse a ′ + b ′ processes and their cross sections σ ab Here, E is the total energy, j a(b) and j ′ a(b) are the initial and final rotational quantum numbers, p 2 a(a ′ )→b(b ′ ) are the initial and final relative momenta between a and b species, σ ab are direct and reverse integral cross sections respectively. The same type of the relationship can be obtained for the thermal rate coefficients k j1j2→j ′ 1 j ′ 2 (T ). Let us rewrite Eq. (8) in the terms of the total energy E. Considering that the relative kinetic energy ǫ between a and b is Eq. (7), we compute the total energy from the lowest possible level between two channels. If it is associated with the second channel a ′ + b ′ , i.e. j ′ 1 j ′ 2 pair, the formula (8) becomes: Here, ] is the energy gap between the direct and reverse channels, and ǫ = p 2 ab /M 12 is the kinetic energy. Now, for the reverse channel the thermal rate coefficient is: Comparing Eqs. (32) and (33) and taking into account Eq. (31) we obtain the detailed balance formula for the thermal rate coefficients: The ratio R j1j2⇋j ′ is proportional to an exponent with the argument depending on ∆E. In this work we computed two direct-reverse processes in the HD + p-H 2 collisions, specifically: with ∆E = 96.6 cm −1 . It would be useful to check how well the computed thermal rates (8) in this work satisfy the detailed balance equation (34). Fig. 9 represents these results. It is shown that all results are in a satisfactory agreement with each other.

IV. CONCLUSION
In this work a rotational method has been applied for the transformation of the 4-dimensional rigid monomer model H 2 -H 2 PES of Ref. 1 to the non-symmetrical potential appropriate for calculations of the HD+H 2 collisions. Different low energy elastic and state-selected inelastic cross sections as well as the thermal rate coefficients for HD+H 2 have been computed and compared with previous calculations, where available. The rotational energy transfer in HD+H 2 is of importance for the thermodynamics of the ISM 24 . By now few and rather conflicting results are available for the low energy HD+H 2 rotational energy transfer, see for example 38 and references therein. The BMKP PES 31 has been already applied to HD+H 2 37,38 . However, this PES may have failures. The fact was mentioned in Refs. 12 and 13 and in 34 . Therefore in this paper a new attempt has been undertaken to carry out alternative computational methods for HD+H 2 collision. In case of the BMKP and DJ PESs the necessary steps for each potential have been described in Secs. II B and II C and also in Ref. 37. In the case of the BMKP potential, which is a full sixdimensional surface 31 , the transformation from H 2 -H 2 to the H 2 -HD system was done by shifting the center of mass in one H 2 molecule to the center of mass of the HD molecule. Because the DJ PES has been formulated for the rigid monomer rotor model, the transformation methodology was more complicated. Simply, the R 1 and R 2 coordinates are not available in this case, they have fixed lengths. In this paper the transformation has been accomplished by rotation of the space fixed OXY Z coordinate system, i.e. by the redirection of the R 3 vector. The new vector R ′ 3 connects the center of masses of the H 2 and the HD molecules as shown in the Fig.  2. This procedure obtains new coordinate angles for the Jacobi vectors R 1 , R 2 , and R 3 , and a new PES as in Eq. (20). New experiments that measure the state-to-state rotational cross sections in the HD+o-/p-H 2 collisions at low temperatures are needed. Thereafter theoreticians and astrophysicists would be able to compare computational results with available experimental data. This type of work was recently accomplished, for instance, for the para-H 2 +H 2 collision 33 . Here it would be useful to mention other contributions on hydrogen-hydrogen collisions [48][49][50] . Additionally, with the use of new HD+H 2 results for the thermal rate coefficients one could carry out new computation of the HD-cooling function mentioned in the introduction 24 .
In conclusion, another interesting system worth mentioning is HD+HD. For this collision there are relatively old experimental state-to-state rotational probabilities for a few selected states 51 . These old data can be useful in comparisons with the computational results obtained with different H 4 PESs: such as the available DJ and the BMKP PESs or some relatively new potentials, for example from works 14,52 . In the case of the DJ surface it would be possible to again apply the rotation procedure of the OXY Z coordinate system as performed in this paper.