© The Authors 2024

1 Introduction

Water, often referred to as the essential source of life, is the most abundant molecule within biological systems and its probable extraterrestrial origin adds an intriguing dimension to its significance (Horneck & Baumstark-Khan 2012). Strong (masertype) emissions from both gaseous and solid forms of water have been observed in many astrophysical environments such as planets, comets, interstellar clouds, and star-forming regions (Ehrenfreund et al. 2003). Interstellar water vapor was discovered in 1969 in the Orion nebula by Charles Townes (Cheung et al. 1969). The discovery of water within the circumstellar disk of the V883 Orionis protostar, observed by the Atacama Large Millimeter/submillimeter Array (ALMA), marks the first instance of water being inherited into a protoplanetary disk without significant changes to its composition (Tobin et al. 2023). Probing water in astrophysical media requires observations of lines that are thermally excited and that do not display population inversion. Unfortunately, such rotational lines are blocked from ground-based observations due to the presence of water in Earth’s atmosphere. Hence, space missions such as the Infrared Space Observatory (ISO), the Submillimeter Wave Astronomy Satellite (SWAS), Odin, Hershel Space Observatory (HSO), and the recently launched mid-infrared space observatory James Webb Space Telescope (JWST) have played a pivotal role in advancing our understanding of water in space. In particular, high-excitation rotational lines of H₂O are commonly detected in the youngest, most deeply embedded protostars, likely excited by far-UV irradiated shocks associated with outflows (Karska et al. 2014). Interpreting these observations requires numerical radiation transfer modeling and relies heavily on accurate collisional rate coefficients. In molecular clouds, the dominant collisional partner for H₂ O is H₂.

The collisional excitation of H₂O by H₂ was the focus of several theoretical studies over the past few decades. The first calculations of rate coefficients for this system were performed by Phillips et al. (1996), considering transitions between rotational levels up to j = 3 in both ortho- and para-H₂O for temperatures ranging from 20 to 140 K.

Subsequent work by Daniel et al. (2011) employed an improved potential energy surface (PES) developed by Valiron et al. (2008), made use of an accurate coupled-channel (CC) method of the MOLSCAT package (Hutson & Sueur 2019), and covered all four symmetries of the H₂O + H₂ system. They computed rate coefficients for transitions involving the first 45 rotational levels (up to j = 11) of ortho-H₂O and para-H₂O, covering a broader temperature range from 5 to 1500 K. For para-H₂O + para-H₂ symmetry, they include rate coefficients for quenching of 44 rotationally excited states of water with rotational levels up to j = 11, collided with hydrogen in its ground rotational state, j = 0. Additionally, for collisions of para-H₂O with para-H₂ in its excited rotational state, j = 2, the rate coefficients are available for the quenching of 19 excited states of water, up to j = 7. For ortho-H₂O + para-H₂ symmetry, they have provided quenching rate coefficients for the excited 44 levels of water collided with both j = 0 and 2 of hydrogen, and additionally for the 9 excited states of ortho-water collided with hydrogen in a j = 4 state. Overall, 8925 individual state-to-state transition rate coefficients are available from the work of Daniel et al. (2011).

In the recent work by Zóltowski et al. (2021), the authors extend the calculations of H₂O + H₂ rotational excitation rate coefficients to temperatures ranging from 10 to 2000 K and provide rate coefficients for 96 rotationally excited states of water with rotational levels up to j = 17, collided with hydrogen treated as a pseudo-atom. They found that to achieve convergence of the cross-sections for H₂O energy levels up to j = 17, it was necessary to include in the rotational basis energy levels up to j = 29. A reduced-dimensionality treatment of the H₂ molecule (that does not fully account for the rotational structure of the H₂ projectile) was adopted to make the scattering calculations computationally feasible. Moreover, to facilitate their calculations even further, they had to appeal to an approximate coupled-states (CS) treatment of scattering, also known as the Coriolis-sudden approximation (Pack 1984; Parker & Pack 1977). Overall, 9312 individual state-to-state transition rate coefficients are available from the work of Zóltowski et al. (2021) Only 1980 of these transitions overlap with the dataset of Daniel et al. (2011).

Furthermore, an important study on H₂O-H₂ collisional rate coefficients was conducted by Faure et al. (2007). They employed quasi-classical trajectory (QCT) calculations to study H₂O-H₂ collisions, including highly-excited states of H₂ and assuming local thermodynamic equilibrium (LTE) populations. Their work provided valuable insights into the behavior of rate coefficients at high temperatures and for rotationally excited H₂. Recent quantum calculations by Stoecklin et al. (2021); Garcia-Vázquez et al. (2024b); García-Vázquez et al. (2024a), Wiesenfeld (2021, 2022), and finally Yang et al. (2024) have investigated the coupling between rotation and bending of H₂O, albeit for a limited set of transitions and kinetic temperatures.

Despite these efforts, there is still a lack of collisional rate coefficients at high temperatures and for highly excited rotational states of H₂O and H₂. Indeed, the only collisional data available for highly excited rotational states of H₂O is that of Zóltowski et al. (2021) computed using an approximate treatment. To address this gap, it is desirable to compute rate coefficients for transitions involving a broader range of collision energies and higher levels of rotational excitation of both collision partners and using alternative methods. In this paper, we present the results of calculations using a rotational basis set that includes 100 states of each para-H₂O and ortho-H₂O (200 states total) and all rotational states of H₂ projectile up to j = 10, for collision energies up to 12000 cm^–1, without invoking the CS approximation; that is, retaining the physics of the Coriolis coupling effect. This became possible due to recent developments of the mixed quantum–classical theory (MQCT) (Bostan et al. 2023, 2024; Joy et al. 2023; Mandal & Babikov 2023; Mandal et al. 2022, 2020) implemented in a user-ready massively parallel computer code MQCT (Mandal et al. 2024; Semenov et al. 2020). Many technical details of these calculations were reported in a recent paper (Joy et al. 2024). Here, we describe how we computed rate coefficients, conduct a rigorous comparison of our results with those available in the literature, and make the dataset of new collisional rate coefficients for the H₂O + H₂ system available to the community.

2 Details of the method

In general, the state-to-state transition rate coefficients are obtained by averaging the corresponding cross sections, ${\sigma }_{n_1{n}_2\to n_1^{\prime }{n}_2^{\prime }}(U)$, over the Boltzmann-Maxwell distribution of kinetic energy, U, at a certain temperature, T. In the MQCT calculations, this procedure has its specifics, which were discussed in detail in a recent paper (Mandal & Babikov 2023). The resultant formula is $$\begin{array}{l}{k}_{n_1{n}_2\to n_1^{\prime }{n}_2^{\prime }}(T)=\frac{v_{\text{ave}}(T)}{{\left(k_{B}T\right)}^2}\frac{e^{-\frac{\mathrm{\Delta }E}{2k_{B}T}}}{\left(2j_1+1\right)\left(2j_2+1\right)}\\ \times {\displaystyle \underset{U_{\min }}{\overset{\infty }{\int }}{\tilde{\sigma }}_{n_1{n}_2{n}_1^{\prime }{n}_2^{\prime }}(U)e^{-\frac{U}{k_{B}T}}\left[1+{\left(\frac{\mathrm{\Delta }E}{4U}\right)}^2\right]\left[1-{\left(\frac{\mathrm{\Delta }E}{4U}\right)}^2\right]U \text{d}U.}\end{array}$$(1)

Here, n₁ n₂ and $n_1^{\prime }{n}_2^{\prime }$ indicate the initial and final states of two colliding molecules. The composite index, $n_1\equiv j_{1k_A{k}_C}$, labels nondegenerate states of the first molecule (using a set of quantum numbers for an asymmetric top, such as water), whereas n₂ ≡ j₂ is used for the second molecule (since this is simply a linear rotor, H₂). The lower limit of integration in MQCT is U_min = |∆E|/4, where $\mathrm{\Delta }E=E_{n_1^{\prime }{n}_2^{\prime }}-E_{n_1{n}_2}$ is the energy difference between the final and initial states of the $n_1{n}_2\to n_1^{\prime }{n}_2^{\prime }$ transition, which is negative for quenching and positive for excitation processes, respectively. v_ave(T) is the average collision speed and k_B is the Boltzmann constant. Finally, ${\tilde{\sigma }}_{n_1{n}_2{n}_1^{\prime }{n}_2^{\prime }}(U)$ is a weighted average of cross sections for excitation and quenching directions of the same transition, introduced to approximately satisfy the principle of microscopic reversibility, as is described elsewhere (Mandal & Babikov 2023): $$\begin{array}{l}{\tilde{\sigma }}_{n_1{n}_2{n}_1^{\prime }{n}_2^{\prime }}(U)=\frac{1}{2}\{\left(2j_1+1\right)\left(2j_2+1\right){\sigma }_{n_1{n}_2\to n_1^{\prime }{n}_2^{\prime }}(U)\\ +\left(2j_1^{\prime }+1\right)\left(2j_2^{\prime }+1\right){\sigma }_{n_1^{\prime }{n}_2^{\prime }\to n_1{n}_2}(U)\}.\end{array}$$(2)

The calculations of collision dynamics were conducted using the AT-MQCT version of the theory, in which the propagation of quantum degrees of freedom is decoupled from the propagation of classical degrees of freedom (Mandal & Babikov 2023). This approach is known to give accurate results for H₂O + H₂. An adaptive-step-size predictor method ADAPTOL was employed for time-propagation, with the tolerance parameter, ε = 10^–3, determined by convergence studies. This method uses larger time steps in the asymptotic region of the PES and then reduces the time step down to 10 a.u. in the interaction region. Initial conditions for MQCT trajectories were sampled randomly using a Monte-Carlo approach with the maximum value of impact parameter, b_max = 25 bohr, for the lowest collision energy and b_max = 15 bohr for all higher energies. These values correspond to the orbital angular momentum quantum numbers l_max = 27 and 285, respectively. Many technical details of these calculations have been reported in a recent paper (Joy et al. 2024).

Our calculations were carried out for ten values of kinetic energy, U = 20.0, 41.28, 84.0, 170.47, 346.41, 703.89, 1430.0, 2906.3, 5906.0, and 12000 cm^–1. For numerical integration, a cubic spline of the entire integrand in Eq. (1) was constructed between these data points for each individual transition and was extrapolated toward the process threshold at low collision energies and toward the high collision energy limit using a smooth analytic function, as is described elsewhere (Mandal & Babikov 2023). The result of this fitting was carefully checked for smoothness before using this dependence for numerical integration over the entire energy range.

Certain astrophysical applications may require the effective rate coefficients, $k_{n_1\to n_1^{\prime }}^{n_2}(T)$, which are obtained by summing state-to-state rate coefficients, $k_{n_1{n}_2\to n_1^{\prime }{n}_2^{\prime }}(T)$, over the final states, $n_2^{\prime }$ of H₂, for a given initial state, n₂ of H₂: $$k_{n_1\to n_1^{\prime }}^{n_2}(T)={\displaystyle \sum _{n_2^{\prime }}{k}_{n_1{n}_2\to n_1^{\prime }{n}_2^{\prime }}}(T).$$(3)

The “thermalized” state-to-state rate coefficient, ${\overline{k}}_{n_1\to n_1^{\prime }}(T)$, between the rotational states of the H₂O molecule can be obtained by averaging over the initial rotational levels, n₂, of para- or ortho-H₂, as is shown below: $${\overline{k}}_{n_1\to n_1^{\prime }}(T)={\displaystyle \sum _{n_2}{w}_{n_2}}{k}_{n_1\to n_1^{\prime }}^{n_2}(T),$$(4)

where the weights of initial states of the projectile H₂ can be expressed as $w_{n_2}=\left(2j_2+1\right)e^{-\frac{E_2}{k_{B}T}}/Q_2(T)$, and the partition function for para- or ortho-hydrogen molecules is $Q_{p/o}\left(T\right)={\displaystyle {\sum }_{n_2}\left(2j_2+1\right)}{e}^{-\frac{E_2}{k_{B}T}}$.

3 Results and discussion

First, we present a comparison of our computed rate coefficients with those from the datasets of Zóltowski et al. (2021) and Daniel et al. (2011) for a subset of 1980 transitions that are present in all three datasets. For our results obtained with the AT-MQCT method (Mandal & Babikov 2023; Mandal et al. 2020), we use the abbreviation “AT.” For the results of Zóltowski et al. (2021), we use the abbreviation “CS,” which corresponds to their CS calculations. For the results of Daniel et al. (2011) we use the abbreviation “CC,” which corresponds to the CC method. We note that all three methods used the same PES of the H₂O + H₂ system (Valiron et al. 2008).

A comparison of three sets of data is conveniently done using the method of analysis called a Dalitz plot (Babikov et al. 2002). For a given transition, the three methods give different rate coefficients, k_AT, k_CS, and k_CC. We defined three unitless variables (coordinates): $$\begin{array}{l}{\zeta }_{\text{AT}}=\frac{k_{\text{AT}}}{k_{\text{AT}}+k_{\text{CS}}+k_{\text{CC}}},\\ {\zeta }_{\text{CS}}=\frac{k_{\text{CS}}}{k_{\text{AT}}+k_{\text{CS}}+k_{\text{CC}}}\text{ and }\\ {\zeta }_{\text{CC}}=\frac{k_{\text{CC}}}{k_{\text{AT}}+k_{\text{CS}}+k_{\text{CC}}},\end{array}$$

which satisfy the relation ζ_AT + ζ_CS + ζ_CC = 1 and which each vary in range from zero to one. These variables were used to place one point (that corresponds to this transition) within the area of a triangular plot that has three axes crossing at 120°. Each point on a Dalitz plot corresponds to one transition and all data points fall within a triangle, as one can see in Figs. 1, 2 and 3. If a data point is close to the middle of the Dalitz plot, this means that we have k_AT ~ k_CS ~ k_CC for this transition (rate coefficients from three dataset are approximately the same). If a point is close to the edge of the triangle, this means that in one of datasets the value of rate coefficient is much smaller than in the other two (e.g., k_CS ≪ k_AT + k_CC). If a point is in the corner of the triangle, this means that in one of datasets the value of rate coefficient is much larger than in the other two (e.g., k_CC ≫ k_AT + k_CS). Finally, if the point is on the diagonal of the triangle and is off-center, it means that only two datasets have similar rate coefficients (e.g., k_AT ≈ k_CC, k_cs).

In Fig. 1, we compare rate coefficients for 990 individual state-to-state transition rate coefficients, $k_{n_1{n}_2\to n_1^{\prime }{n}_2^{\prime }}(T)$, for para-H₂O + para-H₂ collisions at four different temperatures: T = 100 K, 500, 1000, and 1500 K. A similar picture for 990 transitions in ortho-H₂O + para-H₂ is presented in Fig. S1. We note that all these transitions correspond to elastic scattering of the ground state H₂ projectile, (j₂ = 0) → (j₂ = 0), since only these processes are included in the CS dataset of Zóltowski et al. (2021). From these figures, we see that the results of all three methods are in relatively good agreement because the majority of data points appear relatively close to the center of the Dalitz plot. One can notice that at lower temperatures (T = 100 K) the distribution of data points spreads more along the CS axis and remains roughly symmetric around it, but spreads somewhat less along the other two directions. This means that the agreement is a bit better between the AT and CC datasets, and these two sets show about the same differences when compared to the CS dataset. Still, at T = 100 K the center of distribution is in the middle of the plot, meaning that on average all three datasets give similar rates of energy transfer. This picture, however, changes as the temperature rises.

At T = 1500 K, the center of distribution moves off-center, indicating that CS rate coefficients are somewhat higher than those from AT and CC datasets. However, the distribution of data points becomes somewhat much more concentrated, which means that at higher temperatures the agreement improves for most individual state-to-state transitions (with only a few outlying points). Figure S1 shows the same trends for ortho-H₂O.

In Fig. 2, we present the Dalitz plot analysis of the effective rate coefficients, $k_{n_1\to n_1^{\prime }}^{n_2}(T)$, which include summation over the final states of the projectile H₂. For our dataset (AT), the sum over all final states up to j₂ = 8 was computed for all transitions, using Eq. (3). For CC data, the effective rate coefficients from Daniel et al were computed using their code and data retrieved from BASECOL database (Dubernet et al. 2024), which includes final states of H₂ (j₂ = 0, 2, 4). For the CS dataset, where H₂ is not allowed to be excited, state-to-state transition rate coefficients for collisions with H₂ (j₂ = 0) were used instead (i.e., same as in Fig. 1). From Fig. 2, it becomes clear that the effect of rotational excitation of H₂ is quite significant. The majority of data points in this figure, except a few outliers, fall into 1/3 of the Dalitz plot area along one of its edges, indicating a better agreement between the AT and CC datasets where the rotational excitation of H₂ is considered, and a somewhat larger deviation of these two datasets from the CS dataset, where H₂ is described as a pseudo-atom and the value of rate coefficient is considerably smaller. Among the four temperatures presented, this effect is more pronounced at T = 100 K, when the majority of data points are found very close to one of the edges (see Fig. 2), which means that the value of effective rate coefficients in AT and CC dataset are much larger than those in the CS dataset. Figure S2 shows the same trends for ortho-H₂O.

Finally, in Fig. 3 we present the comparison of thermal rate coefficients, ${\overline{k}}_{n_1\to n_1^{\prime }}^{\prime }(T)$, obtained by averaging over the Boltzmann distribution of the initial states of para-H₂ projectiles. For our dataset (AT), the average over even initial states up to j₂ = 8 was computed for all transitions, using Eq. (4). For the CS dataset, where H₂ was not allowed to be excited, state-to-state transition rate coefficients for collisions with H₂(j₂ = 0) were used instead (i.e., same as in Figs. 1 and 2). For the CC dataset, thermal rate coefficients were computed according to the prescription of Daniel et al. (2011) using their data retrieved from the BASECOL database (Dubernet et al. 2024). Namely, for the calculation of thermal rate coefficients Daniel et al. proposed using the effective rate coefficients computed for H₂(j₂ = 2) as an educated guess for H₂(j₂ = 4,6,8) that were not available. This is a very reasonable assumption, which was implemented here for all transitions in ortho-H₂O + para-H₂, as is represented by Fig. S3. For para-H₂O + para-H₂, Daniel et al. (2011) had the H₂(j₂ = 2) data for some transitions, but not for all. For those missing cases, they proposed using the data they had for para- H₂O + ortho-H₂(j₂ = 1) as an approximation for para-H₂O + para-H₂(j₂ = 2, 4, 6, 8), which, again, is a reasonable approximation. We wrote our own code that computes effective rate coefficients using the data published by Daniel et al. (2011) and implements all the assumptions they recommended. This code is available at GitHub¹.

From Fig. 3, we can see that, overall, the thermal rate coefficients of the CS dataset are somewhat smaller than those from the other two datasets and this effect is more pronounced at the lowest temperature considered here, T = 100 K. The agreement between the AT and CC datasets is again somewhat better, and these two sets of data show about the same differences when compared to the CS dataset. Indeed, one can notice that at higher temperatures, T = 1000 and 1500 K, the distribution of points in the Dalitz plot is roughly symmetric around the CS axis. It should be noted, however, that the CS dataset was specifically built for high-temperature applications, and therefore this small disadvantage is not essential. Overall, it appears that, despite approximations made during the calculations of the CS dataset (the CS approximation and the pseudo-atom treatment of H₂), it provides thermal rate coefficients in good agreement with the other two datasets, except perhaps at temperatures of T ∼ 100 K and lower. Figure S3 shows the same trends for ortho-H₂O + para-H₂.

Next, we conducted a detailed comparison of our AT dataset and the CS dataset using 9312 transitions present in both datasets. These include 4656 transitions for quenching of the lower 97 states of para-H₂O. In Figures 4, 5 and 6, we present a comparison of the individual state-to-state $k_{n_1{n}_2\to n_1^{\prime }{n}_2^{\prime }}(T)$, the effective $k_{n_1\to n_1^{\prime }}^{n_2}(T)$, and thermal ${\overline{k}}_{n_1\to n_1^{\prime }}(T)$ rate coefficients, respectively, in the two datasets, at four different temperatures: T = 100, 500, 1000, and 1500 K. Similar pictures for 4656 transitions in ortho-H₂O are presented in Figs. S4, S5 and S6. From Fig. 4 we see, first of all, that a systematically good agreement between the CS dataset (computed using an approximate quantum method) and our AT dataset (computed using an MQCT approach) is found to persist through four orders of magnitude of rate coefficient values, which means that the agreement is good not only for the most intense transitions but for all of them, including those transitions that are ~ 10 000 weaker. At higher temperatures, T = 500, 1000, and 1500 K, the differences between rate coefficients for the individual state-to-state transitions in two datasets (see Fig. 4) are within a factor of two for the majority of transitions. At the lower temperature, T = 100 K, many transitions show somewhat larger differences but, for the group of most intense transitions with rate coefficients larger than 10^–11 cm³ s^–1, the differences between the two datasets remain within a factor of two. Importantly, for this group of transitions, the differences between the two datasets remain within a factor of two at all temperatures. In Table 1, we present the average difference between state-to-state rate coefficients in AT and CS datasets at different temperatures. These values are always within ~50%. For para-H₂O, the largest difference is observed at T = 100 K, with AT rate coefficients being smaller, while for ortho-H₂O the larger difference is observed at higher temperatures, with AT rate coefficients being larger.

From Fig. 5 for the effective rate coefficients, we see that the differences between the two datasets become larger for many transitions, with AT rate coefficients being systematically larger (the same property that we saw in Fig. 2 above), and particularly so at low temperatures. However, in Fig. 5 one can notice that these larger differences mostly concern those transitions that are weaker. For the most intense transitions (with rate coefficients larger than 10^–11 cm³ s^–1), the differences between the effective rate coefficients in the two datasets are usually within a factor of two, with only a few exceptions, and this is the case for all temperatures. The data for ortho-H₂O presented in Figs. S4 and S5 show similar trends.

From Fig. 6 for thermal rate coefficients, we see a very good agreement between the two datasets. Larger differences, present at the level of the effective rate coefficients, disappear at the level of thermal rate coefficients, due to averaging over the distribution of the initial states, which is the same effect we saw in Fig. 3 above. In fact, out of the twelve frames presented in Figs. 4, 5 and 6, the best agreement is found in the last frame of Fig. 4, which corresponds to the comparison of thermal rate coefficients at high temperature, T = 1500 K. Here, for almost all transitions, the differences between rate coefficients in the two datasets are within a factor of two. However, it should be noted that at low temperature, T = 100 K (the leftmost frame of Fig. 6), many transitions show much larger differences and deviate from the main trend that still follows the diagonal of the figure. For these transitions, the rate coefficients of the AT- MQCT dataset are larger than those of the CS dataset. Also, in the leftmost frame of Fig. 6, one can spot a group of points with large values of rate coefficients that exhibit differences exceeding a factor of two. These transitions are: 7₆₂ → 6₅₁, 7₇₁ → 6₆₀, 8₅₃ → 7₄₄, 10₁₉ → 8₁₇, 8₄₄ → 7₃₅, 9₂₈ → 7₂₆, 8₀₈ → 7₁₇, 9₁₉ → 7₁₇, 9₁₉ → 8₀₈, 8₁₇ → 6₁₅, 8₃₅ → 6₃₃, 11_1,11 → 9₁₉. Figure S6 for ortho-H₂O shows similar trends, with only a few transitions (among the intense ones) exhibiting a difference over the factor of two between the two datasets: 7₇₀ → 6₆₁, 11₄₇ → 11₃₈, 8₅₄ → 7₄₃, 8₄₅ → 6₄₃, 8₃₆ → 6₃₄, 9₅₄ → 8₄₅, 7₂₅ → 5₂₃, 9₄₅ → 8₃₆, 6₆₁ → 5₅₀, 11₈₃ → 10₉₂, 8₂₇ → 6₂₅, 7₃₄ → 5₁₄.

Fig. 1 Dalitz plots for the comparison of 990 state-to-state transition rate coefficients for para-H2O + para-H2 collisions at four different temperatures: 100, 500, 1000, and 1500 K. Each triangular plot represents a comparison between three computational methods: MQCT from this work (AT), full quantum (CC) calculations (Daniel et al. 2011), and full quantum (CS) calculations (Zóltowski et al. 2021) using MOLSCAT. Different red dots within each triangle represent different transitions.

Fig. 2 Same as in Fig. 1 but for the effective rate coefficients defined by Eq. (3).

Fig. 3 Same as in Fig. 1 but for the thermal rate coefficients defined by Eq. (4).

Fig. 4 Comparison of 4656 state-to-state transition rate coefficients for quenching of 97 lower states of para-H2O in collision with para-H2 computed using AT-MQCT (this work) vs. those predicted by full quantum CS calculations using MOLSCAT (Zóltowski et al. 2021) at four different temperatures: 100, 500, 1000, and 1500 K. The dashed red lines represent a factor-of-two difference.

4 Implications for astrophysical modeling

While there is a generally good agreement between rate coefficients computed here and those available from the literature, it is important to emphasize the advantage offered by our new dataset, which has an explicit dependence of rate coefficients on the rotational states of H₂ projectile in a broad range. In Fig. 7, we present thermal distribution of the rotational states of H₂ projectile at various temperatures. It demonstrates that an approximation in which the initial states of H₂ are restricted to the ground para-state (j₂ = 0) and ground ortho-state (j₂ = 1) is valid only at ~ 100 K or below, but even in this case it is important to consider ortho-H₂, because its population is significant, about 50% of para-H₂ population (purple curve in Fig. 7). At 500 K, populations of H₂ states j₂ = 2 and 3 become significant (light blue curve in Fig. 7). Importantly, at 1000 K, the population of state j2 = 4 exceeds that of j₂ = 0, and the population of j₂ = 5 is very close to that of j₂ = 0 (green curve in Fig. 7). At 2000 K, populations of all states up to j₂ = 7 exceed that of j₂ = 0 and the population of j₂ = 8 is very close to that (red curve in Fig. 7). Therefore, the effect of these states should not be neglected.

This simple analysis shows that at higher temperatures the rotational states of H₂ up to j₂ = 8 must be taken into consideration. Our dataset is the only one that offers this possibility and computes a true thermal average over the initial states of the projectile H₂ molecules, as is appropriate for temperatures up to 2000 K. Moreover, using the effective rate coefficients, $k_{n_1\to n_1^{\prime }}^{n_2}(T)$, which we provide for various initial states of H₂, users can compute an average over any distribution of the rotational states of the projectile (which may be different from the Boltzmann distribution). This can be necessary in non- LTE conditions, when the distribution of the initial states of H₂ is different from LTE, or when the rotational temperature of the background H₂ gas deviates from the kinetic temperature (Mandal & Babikov 2023). Our dataset offers this unique flexibility as well.

To illustrate these properties, we plot in Figs. 8 and 9 the dependence of effective rate coefficients (for the 20 most intense transitions in para-H₂O) on the initial state of H₂ projectile in the range 0 ≤ j₂ ≤ 9, for several temperatures. Similar data for ortho-H₂O are presented in Fig. S7 and S8. In these figures, for convenience, we plot normalized (or relative) values of the effective rate coefficients, computed as a unitless ratio between the rate coefficient for a certain value of j₂ of H₂ projectile and the same rate coefficient for the ground state of H₂ projectile. Namely, for para-H₂ we plot $R_{j_2}=k_{n_1\to n_1^{\prime }}^{j_2}/k_{n_1\to n_1^{\prime }}^{j_2=0}$at various values of temperature, while for ortho-H₂ we plot $R_{j_2}=k_{n_1\to n_1^{\prime }}^{j_2}/k_{n_1\to n_1^{\prime }}^{j_2=1}$ for the same temperatures. These ratios can be viewed as scaling factors that one would have to apply to the effective rate coefficients computed for the ground state H₂ in order to obtain the effective rate coefficients for the excited states of H₂. Similar moieties were computed in Daniel et al. (2011) for a limited number of transitions with j₂ = 2 and 4. Our analysis is presented for many more transitions and up to j₂ = 9. For para-H₂O, we included the following 20 transitions, identified using Fig. 4 as most intense at 1500 K: 1₁₁ → 0₀₀, 2₁₁ → 2₀₂, 2₂₀ → 1₁₁, 2₂₀ → 2₁₁, 3₁₃ → 2₀₂, 3₁₃ → 2₁₁, 3₂₂ → 3₁₃, 3₃₁ → 2₂₀, 3₃₁ → 3₂₂, 4₀₄ → 3₁₃, 4₁₃ → 4₀₄, 4₂₂ → 4₁₃, 4₃₁ → 4₂₂, 5₁₅ → 4₀₄, 5₂₄ → 4₁₃, 5₃₃ → 4₂₂, 5₃₃ → 5₂₄, 5₄₂ → 4₂₂, and 5₄₂ → 4₃₁. For ortho-H₂O, the 20 most intense transitions are listed in SI.

From Fig. 8, one can see that the largest change in effective rate coefficients happens between the ground and first excited states of para-H₂ (i.e., as we move from j₂ = 0 to 2). This property was observed in the past by several authors who explored the effect of lower excited states of H₂ (Daniel et al. 2011). However, it is often assumed that after that (at j₂ = 4 and above) the rate coefficients remain constant and equal to those with j₂ = 2. Our calculations show that this is not quite true. From Fig. 8, we see that rate coefficients keep growing when we look at j₂ ≥ 4, although ata slower pace. From Fig. 9, we see that a similar trend is valid in the case of ortho-H₂. Namely, the largest increase in effective rate coefficients happens as we move from j₂ = 1 to 3, but the rate coefficients keep growing in the range j₂ ≥ 5 at a slower pace. Similar behavior is observed for ortho-H₂O, as is illustrated by Figs. S7 and S8. On a semi-quantitative level, we conclude that at higher temperatures (1500 K in Figs. 8, 9, S7 and S8), the values of effective rate coefficients for the highest states of H₂ that we considered (j₂ = 8 and 9) are more than twice as large as those for the ground state H₂ (j₂ = 0 and 1), for many transitions. The value of $R_{j_2}=2$ is shown in these figures by a dashed horizontal line, for convenience.

The dependencies of effective rate coefficients on the initial state of H₂ presented in Figs. 8, 9, S7 and S8 are all very similar and are relatively simple, but this simplicity is partially due to the way in which these processes were selected (using Fig. 4 at 1500 K). However, there are different ways to select a small group of representative processes. Therefore, we also tried to identify the 20 most important processes at each individual temperature and for each individual symmetry (of para- or ortho-H₂O collided with para- or ortho-H₂) using the values of effective rate coefficients (like those from Fig. 5, rather than Fig. 4) that include summation over the final states of H₂ and that therefore are more directly related to the final values of thermally averaged rate coefficients. These data are presented in Figs. S9–12. Overall, they demonstrate similar dependencies with rate coefficients growing roughly by a factor of two when the rotational quantum number of H₂ is raised from j2 = 0 and 1 to j2 = 8 and 9, but these dependencies are more complex and show far more variations between different transitions in H₂O, different temperatures, and different symmetries. Therefore, the use of actual data computed in this work, rather than a simplified scaling law that can be easily deduced from Figs. 8, 9, S7 and S8, is expected to give more accurate results, and therefore is recommended.

Fig. 5 Same as in Fig. 4 but for the effective rate coefficients defined by Eq. (3).

Fig. 6 Same as in Fig. 4 but for the thermal rate coefficients defined by Eq. (4).

Fig. 7 Boltzmann distribution of H2 rotational states as a function of rotational quantum number, j2, at five different temperatures between 100 and 2000 K, as is shown by color in the figure. Para-H2 states are represented by filled circles, while ortho-H2 states are shown by empty circles.

Fig. 8 Dependence of effective rate coefficients for para-H2 O on the initial state, j2, of the para-H2 projectile, computed by MQCT for the 20 most intense transitions (identified using Fig. 4) at four temperatures, as is indicated in the figure.

Fig. 9 Same as in Fig. 8 (for the same 20 transitions in para-H2O) but as a dependence on the initial state, j2, of the ortho-H2 projectile.

5 Conclusions

In this paper, we have presented a new dataset of state-to-state, effective, and thermal rate coefficients for transitions between the rotational states of H₂O collided with H₂ background gas. All four symmetries of ortho- and para-water combined with ortho- and para-hydrogen were considered. This dataset offers a significant expansion over the other existing datasets in terms of the rotational states of water (200 states) and the rotational states of hydrogen (10 states). A detailed comparison with previous datasets is presented, which demonstrates that the approximate MQCT employed in this work is sufficiently accurate. The behavior of rate coefficients for collisions with highly excited H₂ molecules is presented for the first time. It shows that collisional rate coefficients for rotational transitions in H₂O molecules grow with the rotational excitation of H₂ projectiles and exceed those of the ground state, H₂, by roughly a factor of two. These findings are important for the accurate description of water molecules in high-temperature environments where the hydrogen molecules of background gas are rotationally excited, and the collision energy is high. The rate coefficients presented here are expected to be accurate up to a temperature of ∼2000 K, when the vibrational excitation of H₂ O bending mode (not considered here) will likely start to play some role. The exclusion of vibrational excitation, particularly the bending mode of H₂O, is a limitation of our study that becomes relevant at higher temperatures. Above 2000 K, we can expect that ro-vibrational coupling could modify rates, new energy transfer pathways involving vibrational excitation or quenching could open up, and the increased internal energy of vibrationally excited H₂O could affect collision dynamics.

The individual state-to-state transition rate coefficients computed in this work will be submitted to the BASECOL database (Dubernet et al. 2024) and the new EMAA (EMAA 2024) database expanding the available resources for astrophysical modeling. The effective rate coefficients are available for a broad range of H₂ rotational states and can be employed for the modeling of collisional energy transfer in nonequilibrium conditions when the distribution of the rotational states of H₂ deviates from the Boltzmann distribution at a given kinetic temperature.

Data availability

Our codes that generate effective and thermal rate coefficients are available from the GitHub site: https://github.com/MarquetteQuantum/MOLRATES.

The supplementary figures are available at Zenodo: https://doi.org/10.5281/zenodo.14171267.

Acknowledgements

This research was supported by NASA, grant number 80NSSC24K0208. D.B. acknowledges the support of the Way Klingler Research Fellowship and the Haberman-Pfletshinger Research Fund. C.J. acknowledges the support of the Schmitt Fellowship. D.B. acknowledges the support of the Eisch Fellowship and the Bournique Memorial Fellowship. We used HPC resources at Marquette funded in part by the National Science Foundation award CNS-1828649.

References

All Tables

All Figures

Fig. 1 Dalitz plots for the comparison of 990 state-to-state transition rate coefficients for para-H2O + para-H2 collisions at four different temperatures: 100, 500, 1000, and 1500 K. Each triangular plot represents a comparison between three computational methods: MQCT from this work (AT), full quantum (CC) calculations (Daniel et al. 2011), and full quantum (CS) calculations (Zóltowski et al. 2021) using MOLSCAT. Different red dots within each triangle represent different transitions.

In the text

	Fig. 2 Same as in Fig. 1 but for the effective rate coefficients defined by Eq. (3).
In the text

	Fig. 3 Same as in Fig. 1 but for the thermal rate coefficients defined by Eq. (4).
In the text

	Fig. 4 Comparison of 4656 state-to-state transition rate coefficients for quenching of 97 lower states of para-H2O in collision with para-H2 computed using AT-MQCT (this work) vs. those predicted by full quantum CS calculations using MOLSCAT (Zóltowski et al. 2021) at four different temperatures: 100, 500, 1000, and 1500 K. The dashed red lines represent a factor-of-two difference.
In the text

	Fig. 5 Same as in Fig. 4 but for the effective rate coefficients defined by Eq. (3).
In the text

	Fig. 6 Same as in Fig. 4 but for the thermal rate coefficients defined by Eq. (4).
In the text

	Fig. 7 Boltzmann distribution of H2 rotational states as a function of rotational quantum number, j2, at five different temperatures between 100 and 2000 K, as is shown by color in the figure. Para-H2 states are represented by filled circles, while ortho-H2 states are shown by empty circles.
In the text

	Fig. 8 Dependence of effective rate coefficients for para-H2 O on the initial state, j2, of the para-H2 projectile, computed by MQCT for the 20 most intense transitions (identified using Fig. 4) at four temperatures, as is indicated in the figure.
In the text

	Fig. 9 Same as in Fig. 8 (for the same 20 transitions in para-H2O) but as a dependence on the initial state, j2, of the ortho-H2 projectile.
In the text