State-to-State Rate Constants for the O(3P)+H2(v) System: Quasiclassical Trajectory Calculations

1. Introduction

Modeling combustion processes in the hydrogen-oxygen system has great practical innovation, since hydrogen is a promising fuel that can improve both the efficiency and environmental friendliness of gas turbine engines and installations. The problem of quantitative description of kinetics in the hydrogen-oxygen system remains unsolved, both with points of agreement between the results of various experiments, and with points of agreement between theoretical models and experiment [1]. The method in presented work to estimate the state-to-state reaction rate constants makes it possible to increase the efficiency of studying hydrogen combustion in advanced gas turbine engines. This approach based on molecular dynamics methods.

Numerical modeling of molecular dynamics is an actively developing field of science, which has become the de facto main method to obtain new information about the rate constants of physicochemical processes. The progress of molecular dynamics methods is due to the availability of powerful computers, the development of computational algorithms (in particular, the increasingly widespread use of neural networks), the ability to perform numerical simulation from the first principles and the absence in a numerical experiment of some natural limitations inherent in laboratory diagnostic methods. In particular, the methods of molecular dynamics allows to obtain the most complete description of the dependence of the probabilities of physicochemical processes on the quantum states (rotational, vibrational, electronic) of reacting particles in a wide temperature range. The complete sets of state-to-state rate constants acquired in this way can be used both to obtain averaged characteristics of the reacting system (for example, thermal rate constants), or used directly in kinetic models for thermally nonequilibrium gas observed, for example, in flows with shock waves [2,3,4,5,6,7,8,9] and detonation waves [10,11,12,13,14], flows in a supersonic nozzle [15,16,17], in gas discharges [18,19,20], during emission/absorption of intense radiation [21,22], in the upper layers of the Earth’s atmosphere [23]. Thus, the development of the methods for state-specific rate constants calculations and the generation of new data on rate constants is an important urgent task, to which this study is devoted.

The main difficulty in estimating the state-to-state rate constants of the reactions and the energy exchange processes is that it is impossible to apply simple estimation approaches (the model of solid spheres, the capture model, the theory of the transition state, etc.) without additional assumptions and empirical corrections. To obtain state-to-state rate constants, it is necessary to use dynamic calculation methods based either on solving the equations of classical mechanics (for example, the method of quasi-classical trajectories (QCT) [24]), or on solving the equations of quantum dynamics of the system under study (for example, the real wave packets method [25]). The latter group of methods is quite difficult to implement even for simple atom-diatom systems, therefore, the main method for estimating the state-to-state rate constants used in most studies is the QCT method [26,27,28,29,30,31,32]. Therefore, in this paper the implementation and application of the QCT method will be consider.

The implementation of the QCT method requires either a preliminary assignment of the entire potential energy surface (PES) of a system of colliding particles in some analytical form (PES can be obtained from ab initio calculations or based on the modelling multidimensional potentials, for example, LEPS potential [33]), or calculating the PES along each calculated trajectory "on the fly" by calculating the energy gradient at the trajectory point (usually by ab initio methods) [34]. In this paper, the first approach was used – the construction of an approximation of the PES based on data from ab initio calculations and subsequent calculations of the trajectory ensembles of interest, since the number of trajectories per ensemble should be large enough so that the relative statistical error in calculating the rate constants is less than 10%.

A separate task is to approximate the PES as a multidimensional function depending on the coordinates of the atomic nuclei of the reacting system, for which approximations based on the least squares method are often used, representing the PES function as a series of degrees of atomic coordinates [35]. An alternative approach is to use a neural network as a universal approximator. One of such methods is approximation by a neural network that accepts symmetrized coordinates as an input (PIP-NN — permutation invariant polynomials-neural network method) [36]. The use of such symmetrized coordinates makes it possible, with the existing indistinguishability in permutations of identical atoms in the reacting system, to reduce the number of PES points necessary for constructing an approximation.

The purpose of the present work is, first, to implement the QCT method for calculating the state-to-state rate constants in the atom-diatom A+BC test system using the PIP-NN method to approximate the PES of such a system, and, second, to obtain a complete set of state-resolved rate constants for the exchange reaction A+BC=AB+C, dissociation and VT (vibrational-translational) relaxation of BC molecules on A atoms. The $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2$ system was chosen as such a test system, since it was previously studied in detail by molecular dynamics methods in the works of other authors [37,38,39,40,41], which will allow to validate our implementation of the QCT method. Nevertheless, in previous studies [37,38,39,40,41], the probabilities of dissociation and VT relaxation of H₂$\left(v\right)$ molecules on the O atom were not considered. So, the QCT method for these processes, as far as the authors know, is used for the first time. The obtained full set of state-to-state rate constants for $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2$ system (including dissociation ant VT relaxation) can be important, for example, for describing the nonequilibrium processes in the Earth’s ozone layer, and approximations for thermal rate constants obtained on its basis can be used in kinetic models of combustion of hydrogen-containing fuels.

2. Methodology

2.1. PES approximation

All calculations of the PES of the $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2$ system along the lower main surface ${}^3{A}^{″}$ were carried out at the level of the extended multi-configuration quasi-degenerate second-order perturbation theory XMCQDPT2 [42]. As an initial approximation of the wave function, the results of the calculation by the multi-configuration self-consistent field method with a dynamic choice of weights DW-CASSCF [43] was used. Such a choice of a multi-configuration calculation method is associated with the need for a smooth construction of the PES, which is most likely provided by the DW-CASSCF method (see Ref. [44] for details). The size of the active space specified in the DW-CASSCF method is 8 electrons distributed in 6 orbitals, i.e. the full valence space for the considered system. The main set of basis functions in all calculations is the family of Dunning basis sets with diffuse functions aug-cc-pVXZ (X=D, T, Q) [45]. Based on calculations in the basis sets aug-cc-pVXZ (X=D, T, Q), the energies were calculated at each point of the PES in the limit of the complete basic set aug-cc-pV∞Z using the Riemann zeta function [46]. All ab initio calculations were carried out in the Firefly QC v.8.2.0 software package [47], partially based on the GAMESS (US) [48] source code.

The obtained set of PES data was approximated using the PIP-NN method [36]. This approach to fit PES has important advantages: the universality of the method, the availability of neural network training tools (for example, TensorFlow [49], PyTorch [50], Keras [51] software packages, etc.), the availability to implement the analytical nuclear gradient at any point of the PES. At the same time, the main drawback of this approach is related to the peculiarity of the neural network device — for a neural network with a sigmoid activation function or a hyperbolic tangent, the main calculation time of the output neuron is spent on multiple calculations of this function, as well as matrix-vector multiplication operations.

In the framework of the PIP-NN method, before constructing and training a neural network that plays the role of a universal approximator, it is necessary to transform the Cartesian coordinates of atomic nuclei into symmetrized coordinates. Since the $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2$ system corresponds to the AB₂ type, the general form of symmetrized coordinates is already known via Morse coordinates [52]:

$$\left\{\begin{array}{c}{G}_1=exp\left(-\lambda r_{{\mathrm{H}}_1{\mathrm{H}}_2}\right)\hfill \\ G_2=exp\left(-\lambda r_{\mathrm{O}{\mathrm{H}}_1}\right)+exp\left(-\lambda r_{\mathrm{O}{\mathrm{H}}_2}\right)\hfill \\ G_3=exp\left(-\lambda r_{\mathrm{O}{\mathrm{H}}_1}\right)·exp\left(-\lambda r_{\mathrm{O}{\mathrm{H}}_2}\right)\hfill \end{array}\right.$$

(1)

Here $\lambda$ is a dimensional coefficient equal to 1-2 Å⁻¹. In general, it is necessary to solve a combinatorial problem with respect to permutations of atoms and the corresponding permutations of interatomic distances. Software packages such as MAGMA [53] and SINGULAR [54], which are able to determine permutation groups based on a given permutation matrix, help in solving these problems. The coordinates obtained in this way within the PIP-NN approach are fed to the input of the neural network, and at the output the neural network gives the potential energy value.

In this paper, to approximate the PES of the $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2$ system, a neural network of the multilayer perceptron type with two hidden layers was trained, 45 neurons being on the first hidden layer, 55 on the second. The RMS error after training was about 10 meV.

2.2. Rate constants estimation

The direct calculation of the rate constants by the QCT method is based on a statistical analysis of the outcomes of trajectory ensembles with specified conditions. We will briefly describe the stages of all calculations that eventually lead to estimates of the rate constants of interest, according to the methodology from [24,55,56].

First of all, it is necessary to convert the Cartesian coordinates of the atomic nuclei into Jacobi coordinates for the atom–diatomic molecule system. In the case of the $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2$ system, these are the three coordinates of the H₁ atom relative to the H₂ atom, the three coordinates of the $\mathrm{O}{(}^{3}P)$ atom relative to the center of mass of the H₁ H₂ molecule, and the three coordinates of the center of mass of the entire $\mathrm{O}{(}^{3}P)$H₁H₂ system. After the corresponding linear transformation, the Hamiltonian of the system will take the form:

$$H\left(\mathbf{p},\mathbf{q}\right)=\sum _{i=1}^3{\displaystyle \frac{p_i^2}{2{\mu }_{{\mathrm{H}}_1{\mathrm{H}}_2}}}+\sum _{i=4}^6{\displaystyle \frac{p_i^2}{2{\mu }_{\mathrm{O},{\mathrm{H}}_1{\mathrm{H}}_2}}}+\sum _{i=7}^9{\displaystyle \frac{p_i^2}{2M_{\mathrm{O}{\mathrm{H}}_1{\mathrm{H}}_2}}}+V\left(\mathbf{q}\right).$$

(2)

Here $\mu$ are appropriate reduced masses, $M_{\mathrm{O}{\mathrm{H}}_1{\mathrm{H}}_2}$ is the total mass of the entire system, $\mathbf{q}$ or $q_i$ are the Jacobi coordinates, and $\mathbf{p}$ or $p_i$ are the generalized momentum corresponding to $\mathbf{q}$, $V\left(\mathbf{q}\right)$ is the potential energy of the system at the point with the coordinates $\mathbf{q}$. The resulting Hamiltonian (2) allows to compose the equations of motion of Hamiltonian mechanics, which were then solved numerically using an explicit one-step Runge-Kutta-Merson method of the 4th order of accuracy with a variable integration step.

Each ensemble of trajectories was set by initial conditions, expressions for which are described in detail in [24,55,56]. Here we only note that the calculation of the reaction cross section was carried out without scanning over the impact parameter, but each ensemble of trajectories was set only by the vibrational and rotational numbers of the H₂ molecule and the kinetic energy of the incoming $\mathrm{O}{(}^{3}P)$ atom.

Also, to construct the initial conditions and for subsequent analysis of each trajectory, it is necessary to know the vibrational-rotational levels of H₂ and OH molecules, but the approximation of the anharmonic Morse potential may be inaccurate when describing upper vibrational-rotational levels. Therefore, the energy levels of diatomic molecules were determined by numerically solving the radial component of the stationary Schrödinger equation for a diatomic molecule with the LEVEL16 program [57]. However, it is worth noting that this program, alongside with bound states, also calculates the energy levels of quasi-bound states. Despite the fact that some of the quasi-bound states may be metastable and have a lifetime of more than a day, such states will not be considered within the framework of the QCT method, and the trajectories that led to such an outcome will be attributed to the process of dissociation.

Thus, it becomes possible to calculate the cross sections of the following reactions in this system:

$$\begin{array}{c}\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2(v_1,J_1)\to \mathrm{H}+\mathrm{OH}(v_2,J_2),\\ \mathrm{O}{(}^{3}P)+{\mathrm{H}}_2(v_1,J_1)\to \mathrm{O}{(}^{3}P)+{\mathrm{H}}_2(v_2,J_2),\\ \mathrm{O}{(}^{3}P)+{\mathrm{H}}_2(v,J)\to \mathrm{O}{(}^{3}P)+\mathrm{H}+\mathrm{H}.\end{array}$$

(3)

The reaction cross section was calculated using the ratio of the number of trajectories $N_{\mathrm{R}}$ that led to the required outcome to the total number of trajectories of the ensemble N according to the formula from [24,55,56]:

$$\sigma \left(v,J,E\right)\approx \pi b_{\max }^2{\displaystyle \frac{N_{\mathrm{R}}}{N}}$$

(4)

with the statistical error

$$\Delta \left(v,J,E\right)\approx \pi b_{\max }^2{\displaystyle \frac{N_{\mathrm{R}}}{N}}\sqrt{{\displaystyle \frac{N-N_{\mathrm{R}}}{N{N}_{\mathrm{R}}}}}.$$

(5)

Here, $b_{\max }$ is the maximum value of the impact parameter at which the probability of a particular process (reaction, dissociation or VT relaxation process) is above zero, E is the kinetic energy of the incoming $\mathrm{O}{(}^{3}P)$ atom. In all calculations, the total number of trajectories that compose one ensemble was equal to 20000, and the highest kinetic energy of the $\mathrm{O}{(}^{3}P)$ atom was equal to 9 eV.

On the basis of the obtained cross sections the rate constants $k(v,J,T)$ and their error $\Delta k(v,J,T)$ were calculated by the formulas:

$$\begin{array}{c}k(v,J,T)=N_{\mathrm{A}}{\displaystyle \frac{g_{\mathrm{PES}}^{\mathrm{e}}}{g_{\mathrm{R}}^{\mathrm{e}}}}\sqrt{{\displaystyle \frac{8k_{\mathrm{B}}T}{\pi {\mu }_{\mathrm{A},\mathrm{BC}}}}}\underset{0}{\overset{\infty }{\int }}\sigma \left(v,J,E\right)exp\left(-{\displaystyle \frac{E}{k_{\mathrm{B}}T}}\right){\displaystyle \frac{dE}{{\left(k_{\mathrm{B}}T\right)}^2}},\\ \Delta k(v,J,T)=N_{\mathrm{A}}{\displaystyle \frac{g_{\mathrm{PES}}^{\mathrm{e}}}{g_{\mathrm{R}}^{\mathrm{e}}}}\sqrt{{\displaystyle \frac{8k_{\mathrm{B}}T}{\pi {\mu }_{\mathrm{A},\mathrm{BC}}}}}\underset{0}{\overset{\infty }{\int }}\Delta \left(v,J,E\right)exp\left(-{\displaystyle \frac{E}{k_{\mathrm{B}}T}}\right){\displaystyle \frac{dE}{{\left(k_{\mathrm{B}}T\right)}^2}}.\end{array}$$

(6)

Here $g_{\mathrm{PES}}^{\mathrm{e}}$ and $g_{\mathrm{R}}^{\mathrm{e}}$ are the electronic statistical weights of the studied PES and reagents, respectively, $k_{\mathrm{B}}$ is the Boltzmann constant. The integrals of formulas (6) were calculated numerically with piecewise linear interpolation of the reaction cross section $\sigma \left(v,J,E\right)$:

$$\sigma \left(v,J,E\right)=a_{i}E+b_i,$$

(7)

where $a_i={\displaystyle \frac{{\sigma }_i-{\sigma }_{i-1}}{E_i-E_{i-1}}}$ and $b_i={\displaystyle \frac{{\sigma }_{i-1}{E}_i-{\sigma }_i{E}_{i-1}}{E_i-E_{i-1}}}$. Such an approximation of the reaction cross section makes it possible to calculate the rate constant $k\left(T\right)$ by simple summation according to the following formula (the error of such a numerical approach is many times less than the statistical error of trajectory calculations, usually amounting to 5-15% [55]):

$$\begin{array}{c}k(v,J,T)=N_{\mathrm{A}}{\displaystyle \frac{g_{\mathrm{PES}}^{\mathrm{e}}}{g_{\mathrm{R}}^{\mathrm{e}}}}\sqrt{{\displaystyle \frac{8k_{\mathrm{B}}T}{\pi {\mu }_{\mathrm{A},\mathrm{BC}}}}}{\displaystyle \underset{0}{\overset{\infty }{\int }}}\sigma \left(v,J,E\right)exp\left(-{\displaystyle \frac{E}{k_{\mathrm{B}}T}}\right){\displaystyle \frac{dE}{{\left(k_{\mathrm{B}}T\right)}^2}}\approx \\ N_{\mathrm{A}}{\displaystyle \frac{g_{\mathrm{PES}}^{\mathrm{e}}}{g_{\mathrm{R}}^{\mathrm{e}}}}\sqrt{{\displaystyle \frac{8}{\pi {\mu }_{\mathrm{A},\mathrm{BC}}{k}_{\mathrm{B}}T}}}{\displaystyle \sum _{i=0}^n}\left[exp\left(-{\displaystyle \frac{E}{k_{\mathrm{B}}T}}\right)\left(a_i\left(E^2+2E{k}_{\mathrm{B}}T+2{\left(k_{\mathrm{B}}T\right)}^2\right)+b_i\left(E+k_{\mathrm{B}}T\right)\right)\right]{|}_{E_{i-1}}^{E_i}.\end{array}$$

(8)

From the obtained vibrational-rotational state-to-state rate constants (6) it is possible to determine the vibrational state-resolved rate constants $k_v\left(T\right)$ (assuming that the translational and rotational degrees of freedom of the reagents and the products are in equilibrium) and the thermally equilibrium rate constant $k_{\mathrm{term}}\left(T\right)$ as follows:

$$\begin{array}{c}{k}_v\left(T\right)={\displaystyle \sum _{J=0}^{M\left(v\right)}}w(v,J,T)k(v,J,T)\\ k_{\mathrm{term}}\left(T\right)={\displaystyle \sum _{v=0}^N}{\displaystyle \sum _{J=0}^{M\left(v\right)}}w(v,J,T)k(v,J,T),\end{array}$$

(9)

where $w(v,J,T)$ is the normalized population of the H₂ energy level with the vibrational and rotational numbers v and J, N is the number of vibrational levels of the molecule H₂, $M\left(J\right)$ is number of rotational levels in the v-th vibrational level.

3. Results and discussion

The methodology described above for obtaining state-to-state rate constants and its specific implementation naturally needs validation. As validation data, the results of QCT calculations of rate constants or calculations using direct modeling of quantum dynamics performed by other authors, as well as data on rate constants borrowed from experimental works can be used. For the considered $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2$ system, the most widely presented in the literature are the state-to-state rate constants for the $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2(v,J)\to \mathrm{H}+\mathrm{OH}$ process. So, in Figure 1 the cross sections for the $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2(v,J=0)\to \mathrm{H}+\mathrm{OH}$ reaction from our QCT estimates are presented in comparison with the similar QCT results from [38]. One can notice a slight overestimation of cross sections in our modelling compared to [38]. This may be due to the fact that the energy barrier of the reaction $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2\to$ H+OH within the framework of the transition state theory for the PES obtained in this work turned out to be slightly lower ($E_{\mathrm{a}}$=11.5 kcal/mole) than the corresponding value from [38] ($E_{\mathrm{a}}$=13 kcal/mole).

Figure 2 shows a comparison of the state-to-state rate constants found by the QCT method in this work and obtained on the basis of quantum dynamics calculations in [41]. The above comparison shows a good qualitative agreement of the results of the two computational studies. Nevertheless, for reaction channels with a quantum number $v>4$, the difference in the values of the rate constants determined by different methods can be significant, and in the temperature range $T<500$ K can reach an order of magnitude. This discrepancy may be due to the fact that the QCT method cannot take into account quantum effects, such as the tunnel effect and the over-barrier reflection.

A comparison of the thermally equilibrium $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2\to \mathrm{H}+\mathrm{OH}$ reaction rate constant, calculated in this work via the QCT method and obtained in [41] on the basis of quantum dynamics calculations, is shown in Figure 3. In addition, Figure 3 shows an estimate of this rate constant obtained using the non-variational transition state theory adjusted for the tunnel effect and various Arrhenius approximations of the rate constant used in the kinetic models [58,59,60]. One can notice an expected overestimation of the rate constant by the transition state theory, while the QCT method gives a closer values, although somewhat higher at $T<2500$ K, relative to the known data. Such an overestimation has the same cause as in the case of the reaction cross section — lowering the energy barrier at the XMCQDPT2/aug-cc-pV∞Z level of theory. The discrepancy in the values of the rate constant at $T>2500$ K is explained by the different temperature ranges for which the Arrhenius approximation was obtained. In the temperature range $T=100-5000$ K, the thermally equilibrium rate constant of the exchange reaction obtained in this work can be approximated with good accuracy by the Arrhenius dependence (cm³/(mole s)):

$$k_{\mathrm{exch}}\left(T\right)=1.21·{10}^9{T}^{1.37}exp\left(-3684/T\right).$$

(10)

For the dissociation reaction and the VT relaxation process of H₂ molecules in collision with O atoms, there is no data on rate constants in the literature (as far as the authors know). Nevertheless, in order to validate the calculation method, it is interesting to compare the estimates of the rate constants of these processes obtained by the QCT method with the available information on the rates of dissociation and VT relaxation of hydrogen molecules in collisions with other partners. For the dissociation process, such a comparison is presented in Figure 4. It turns out that the equilibrium rate constant of the dissociation of H₂ obtained in this work significantly exceeds those with the different partners available in the literature. However, at temperatures $T>1000$ K, the dissociation rate constant of hydrogen upon collision with O atoms coincides by the order of magnitude with the dissociation rate constants in collision with H and H₂O. The thermally equilibrium rate constant of H₂ dissociation obtained in this work in the temperature range $T=100-5000$ K can be described with good accuracy using the double Arrhenius dependence (cm³/(mole s)):

$$k_{\mathrm{diss}}\left(T\right)=1.92·{10}^9{T}^{0.63}exp\left(-39143/T\right)+1.21·{10}^{15}{T}^{0.19}exp\left(-47567/T\right).$$

(11)

Figure 5 shows a comparison of the VT relaxation time VT of H₂ molecules in collisions with various partners. Note that for the processes of hydrogen relaxation on H, H₂ and O, the relaxation time was determined using the VT rate constant of the transition from the first excited vibrational level to the ground state:

$$p\tau ={\displaystyle \frac{k_{\mathrm{B}}T}{k_{\mathrm{V},\mathrm{T}}^{1,0}\left(1-exp\left(-{\theta }_v/T\right)\right)}}.$$

(12)

Here $k_{\mathrm{V},\mathrm{T}}^{1,0}$ is the VT relaxation rate constant, ${\theta }_v$ is the characteristic vibrational temperature. It can be seen from Figure 5 that according to the results of calculations performed in this work, the vibrational relaxation of H₂ proceeds faster on the O atom than on the Ar atom, but slower than on the H atoms, what is consistent with the qualitative expectations of the dependence of the relaxation time of the molecule on the mass of the collisional partner, expressed, for example, in the well-known Millikan–White formula [65]. At the same time, the results of the calculation using the QCT method predict a more effective VT relaxation of H₂ in collisions with O atoms than in collisions with molecular hydrogen. It can be assumed that this effect is due to the fact that a molecule as a collisional partner, in addition to translational degree of freedom, has rotational and vibrational ones, unlike a single atom.

The lack of data in the literature on the probabilities of dissociation and VT relaxation of H₂ molecules in collision with O atoms is due to the fact that the exchange reaction O+H₂$\to$ H+OH is considered as a more probable outcome of the collision of H₂ and O and hence is more important for kinetic models. The rate constants calculated in this work for these processes are shown in Figure 6. The comparison shows that for H₂ molecules in the ground vibrational state, the exchange reaction is indeed a much more probable process than dissociation. However, for vibrationally excited H₂ molecules, the rate constants of dissociation and VT relaxation upon collision with O can be comparable or even many times higher than the corresponding exchange reaction constants. Thus, the data obtained in this paper on the state-to-state rate constants for dissociation and VT relaxation of molecular hydrogen can be an important addition for the kinetic models describing hydrogen oxidation at high translational temperatures or under essentially nonequilibrium conditions.

4. Conclusion

In this study, an implementation of the QCT method was created for atom–diatom systems, the feature of which is the approximation of the PES of the system by a neural network. Using this implementation, cross sections and state-to-state rate constants for exchange reactions, dissociation and VT relaxation processes were calculated for the $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2$ system. Based on the data on the state-to-state rate constants for exchange and dissociation reactions, Arrhenius approximations of thermally equilibrium rate constants over a wide temperature range were obtained.

It was shown that the cross sections and rate constants obtained in this work for the exchange reactions $\mathrm{O}{(}^{3}P)+{\mathrm{H}}_2(v,J=0)\to \mathrm{H}+\mathrm{OH}$ are in a good qualitative agreement with the results of the QCT and the quantum dynamics method calculations presented in the works of other authors [37,38,41]. In addition, the thermally equilibrium rate constant obtained for the exchange reaction is in good agreement with the corresponding rate constants used in various kinetic models [58,59,60].

As far as the authors are aware, there are no reliable data on the rate constants of dissociation and VT relaxation of H₂ molecules in collision with O atoms in the literature, and the study of these processes by the QCT method was performed in this work for the first time. A comparative analysis of the obtained state-to-state rate constants showed that when H₂ on the lower vibrational levels ($v<5$) collides with an O atom, the probability of its dissociation is negligible compared to the probability of an exchange reaction (however, this difference decreases with increasing temperature). But for the upper vibrational levels of hydrogen molecule, the rate constants of dissociation and VT relaxation of H₂ on O are comparable or may even exceed the corresponding state-resolved rate constants of the exchange reaction. Thus, the rate constants of dissociation and VT relaxation of H₂ on O obtained in this work can be useful in constructing kinetic models describing hydrogen oxidation under thermally nonequilibrium conditions or at high temperatures.

The results of this work make it possible to increase the efficiency of modeling processes in gas-turbine power plant using promising types of fuel, such as methane-hydrogen mixtures or pure hydrogen, including at the high combustion temperature inherent in hydrogen fuels.

In the future, it is planned to refine the calculation of the rate constants of energy exchange processes in this system, since in this work the interaction of the $\mathrm{O}{(}^{3}P)$ atom and the H₂ molecule along the other two terms of this system was not taken into account. In addition, due to the versatility of our PES approximation method, the created computational complex can be easily adapted for other triatomic systems and other configurations of collisions, for example, for a recombination reaction A+B+C.