Introduction

Inelastic collisions between water molecules, H 2 O + H 2 O, play an important role in collisional energy transfer in environments such as atmospheres of icy planets (Wirström et al. 2020;Vorburger et al. 2022), in cometary comae (Cochran et al. 2015;Cordiner et al. 2022), and in the atmospheres of waterrich exoplanets (García Muñoz et al. 2024). In nonequilibrium conditions, which often exist in these environments, the collisional energy transfer competes with radiative energy transfer, alters the populations of molecular states, and affects the spectral properties of the emitted radiation (submillimeter, terahertz, and deep infrared). Therefore, an interpretation of the radiation that is emitted or absorbed by molecules in the environments that deviate from local thermodynamic equilibrium (not in local thermodynamic equilibrium, NLTE) requires detailed modeling using codes such as RADEX (Van der Tak et al. 2007), LIME (Brinch & Hogerheijde 2010), or MOLPOP (Ramos & Elitzur 2018), which in turn require as input accurate rate coefficients for collision-induced state-to-state transitions in the molecules of interest. For this purpose, the databases of collisional rate coefficients such as BASECOL (Dubernet et al. 2024) and LAMDA (van der Tak et al. 2020) were developed, and active computational research is ongoing to populate these databases with more rate coefficients to include new molecular collision partners and to broaden the temperature range. The prediction of the rate coefficients requires a quantum mechanical treatment of the internal molecular states, which is implemented in computational chemistry programs such as MOLSCAT (Hutson & Le Sueur 2020, 2019) and HIBRIDON (Alexander et al. 2023), and more recently, MQCT (Mandal et al. 2024).

The H 2 O + H 2 O collision process is very difficult to describe with quantum mechanics (Agg & Clary 1991;Boursier et al. 2020;Buffa et al. 2000) because both collision partners must be treated as general asymmetric top rotors, which creates a very large number of individual quantum states of the H 2 O + H 2 O system as a whole. Recently, noticeable progress in the modeling of H 2 O + H 2 O rotationally inelastic collisions was achieved using a mixed quantum/classical theory approach (Mandal et al. 2023) that is implemented in the code MQCT. This approximate method treats the relative translational motion of collision partners classically (using the mean-field trajectories), while the internal motion of collision partners such as rotations and vibrations are treated with quantum mechanics (time-dependent Schrödinger equation). Based on this approach, first, a database of thermally averaged cross sections (TACS) was reported (Mandal & Babikov 2023a), which was then converted into a database of thermal rate coefficients (Mandal & Babikov 2023b) for 231 transitions in para-H 2 O and 210 transitions in ortho-H 2 O due to collisions with a thermalized bath of H 2 O molecules in the temperature range 5 ≤ T ≤ 1000 K.

It should be stressed, however, that all databases that were developed so far for the H 2 O + H 2 O system involve averaging over the rotational states of the projectile H 2 O assuming a Boltzman distribution of these states at given kinetic temperature, that is, the equilibrium conditions. These data can probably be used for NLTE simulations, but only when the deviation from LTE is relatively small. Therefore, it is still desirable to have a database of state-to-state transitions in H 2 O + H 2 O that is specifically built for NLTE simulations.

In principle, a database of rate coefficients might be computed at a certain kinetic temperature T kin for all possible state-to-state transitions in H 2 O + H 2 O system, including all individual states of the target and projectile molecules (i.e., without averaging over the distribution of the internal states of the projectile H 2 O). Then, in the modeling of radiation transfer, these rate coefficients might be averaged over the states of the projectile H 2 O using a non-Boltzmann distribution of the internal states (with the appropriate distribution for the NLTE conditions of interest). The practical realization of this approach is nearly impossible for two reasons, however. First, the number of individual transitions in H 2 O + H 2 O system is unmanageably large, close to 1.3 million. Second, even if such a database were prepared, it could not be used because the standard radiative transfer codes mentioned above are not equipped to take these state-specific rate coefficients as input for the modeling.

An alternative approach, which is undertaken in this paper, is to assume that although there is no equilibrium between all possible degrees of freedom in the system, the translational and rotational degrees of freedom can be expected to achieve their own equilibria independently, and that they can be approximately characterized by Boltzmann distributions at two different temperatures, T kin and T rot . Then, T kin can be used to average over the collision energies of two H 2 O molecules, while T rot can be used to average over the rotational states of the projectile H 2 O. This approximation permits us to reduce the dataset to a manageable size, but also it causes the rate coefficients for state-to-state transitions in the target H 2 O to become dependent on two variables: T kin and T rot (instead of one common temperature in the case of equilibrium). This method was applied to CO + CO in cometary coma and to H 3 O + , H 2 O, and HDO collided with H 2 in the interstellar medium (Cordiner et al. 2022;Demes et al. 2023;Faure et al. 2024). A similar two-temperature model is widely used for the description of electrons and ions in nonequilibrium plasma of a gas discharge (Raizer & Allen 1997) because the collisional energy exchange between light and heavy particles is inefficient.

In this paper, we outline a method that permits us to implement this idea in the case of H 2 O + H 2 O, and we explore the implications of nonequilibrium conditions for an astrophysical modeling of H 2 O-rich environments, such as the atmospheres of icy planets and cometary comae.

Details of the method

The standard procedure for predicting thermal state-to-state transition rate coefficients k n 1 →n ′ 1 (T ) begins with averaging the state-to-state transition cross sections

(1)

In this expression, v ave (T ) = 8k B T/πµ is the average collision velocity, k B is the Boltzmann constant, µ is the reduced mass of the two collision partners, and subscripts n 1 n 2 and n ′ 1 n ′ 2 indicate their initial and final states, respectively. Each n is a composite index that represents a full set of quantum numbers for one molecule. Namely, for rotational states of an asymmetric-toprotor molecule, such as water, n denotes j k A k C . If we focus on the first water molecule and treat it as a target, while the second water molecule is treated as a projectile, the next step is the summation over the final states and averaging over the initial states of the projectile,

where the thermal populations (weights) of the initial states of the projectile are defined as

Here E 2 denotes the energies of the rotational states n 2 of the projectile, and Q 2 (T ) is the standard partition function of the projectile,

Since there in no ortho-to-para conversion (of the projectile states n 2 ) during nonreactive inelastic collisions (Lique et al. 2014), Eqs. ( 2)-( 4) are applied separately to the para-and orthostates of the projectile as if they would be different species, which gives two values of k n 1 →n ′ 1 (T ) for each n 1 → n ′ 1 transition in the target molecule. These need to be added together considering the nuclear spin weights of the ortho-and para-states of the projectile molecule. For example, for H 2 O water molecules, the ortho-to-para ratio is 3 to 1, and the weights are 1/4 and 3/4 for contributions to k n 1 →n ′ 1 (T ) of the para-and ortho-states of the projectile.

In nonequilibrium conditions, this procedure is modified as follows. It is assumed that the translational (relative) motion of the molecules is thermalized at certain kinetic temperature T kin , and this temperature enters the analog of Eq. ( 1), where the averaging over molecule-molecule collision energies is computed,

Then, it is assumed that although the rotational motion is also thermalized, it is characterized by a different temperature T rot , and T rot T kin . This rotational temperature describes the distribution of the internal states of the projectile,

and is used for averaging over this distribution,

In this way, the value of state-to-state transition rate coefficient depends on two temperatures, T kin and T rot . At equilibrium, when A208, page 2 of 7

T rot = T kin , the result of Eqs. ( 5)-( 8) is equivalent to that of Eqs. ( 1)-( 4), which can be used to confirm the correctness of the code and the resultant rate coefficients.

The practical implementation of this approach encounters methodological difficulties in the case of a complex molecular system, such as H 2 O + H 2 O. The reason is the huge number of individual state-to-state transitions n 1 n 2 → n ′ 1 n ′ 2 in a model that considers the rotational states of both target and projectile molecules. For example, as previously reported (Mandal & Babikov 2023a), we focused on 231 para-para and 210 orthoortho lower-energy transitions in the target H 2 O, and we used a rotational basis set of 38 para-and 38 ortho-states for the projectile H 2 O molecule (2888 transitions in the projectile, excluding ortho-to-para conversion), which led to close to 1.3 million transitions overall. For each of these, the energy dependence of the cross section σ n 1 n 2 →n ′ 1 n ′ 2 (E), computed on a grid of points, needs to be smoothly interpolated between these points and extrapolated toward the process threshold at low collision energies and toward the high-energy limit using a smooth function. The result of this fitting must be confirmed in some way (e.g., visually) before this dependence is numerically integrated over the collision energy E. While this is possible for a small number of transitions in simple systems, this procedure becomes unfeasible when the number of transitions is about one million or larger, as is the case here.

Therefore, we developed a different approach. We substituted Eq. ( 5) into Eq. ( 8) and exchanged the order of integration over collision energies (taking it outside) with a summation over the states of the projectile (bringing it inside), which gives

We introduced a thermally averaged cross section (TACS) for the transition n 1 → n ′ 1 , defined as

The TACS are computed from the original state-to-state transition cross sections as the sum over the final and average over the initial states of the projectile at chosen values of rotational temperature T rot and for all values of the considered collision energy E, making each TACS a function of these two variables. The number of transitions between the states of a target molecule is relatively small. For example, in our case, it was 231 transitions in para-water and 210 transitions in ortho-water (for the quenching and the corresponding excitation processes), which were computed at six values of the collision energy. Therefore, it is quite manageable to check the behavior of all TACS before they are integrated in Eq. ( 9). Equation ( 9) can be used as is to process the results of exact methods, such as full quantum scattering theory (Demes et al. 2023). We built upon the method developed by Mandal & Babikov (2023b) specifically for the mixed quantum/classical calculations, which permitted us to enforce microscopic reversibility using the cross-section calculations in the excitation and quenching directions. In this case, Eq. ( 9) is replaced by a slightly more complicated formula 11) for several transitions between the ground state 0 00 of para-H 2 O and its several low-lying rotational states (as indicated by different colors in the figure) at T kin = 100 K. The solid, dashed, and dash-dotted lines correspond to a rotational temperature of T rot =10, 50, and 100 K, respectively. (Mandal & Babikov 2023b),

Here, U is the kinetic energy of the mixed quantum/classical trajectories, and U min = ∆E/4 is its threshold value for a given n 1 → n ′ 1 transition. ∆E = E n ′ 1 -E n 1 is the energy change in the transition, which is negative for quenching and positive for excitation processes. Finally, the TACS in Eq. ( 11) represents a weighted average of the TACS for the quenching and excitation directions of the same transition (Mandal & Babikov 2023b),

For a detailed derivation of Eqs. ( 11)-( 12), we refer to Mandal & Babikov (2023b). Each TACS in Eq. ( 12) was computed as the sum over the final and average over the initial states of the projectile, as before,

Results and discussion

From previous work (Mandal & Babikov 2023a), the crosssections σ n 1 n 2 →n ′ 10 -12 10 -11 10 -10 10 -9 10 -8 10 -12 10 -11 10 -10 10 -9 10 -8 k (cm 3 s -1 ) T = 200 K 10 -12 10 -11 10 -10 10 -9 10 -8 k (cm 3 s -1 ) T = 500 K 10 -12 10 -11 10 -10 10 -9 10 -8 T = 800 K Fig. 3. Same as Fig. 2, but for the quenching of 210 states of ortho-H 2 O.

T rot using Eq. ( 10). The excitation and quenching cross sections were combined in Eq. ( 12) to give σn 1 n ′ 1 (U, T rot ). These were interpolated and integrated over U for any given input temperature T kin according to Eq. ( 11) using a numerical procedure (Mandal & Babikov 2023b) with minor modifications. Figure 1 gives examples of the function under the integral of Eq. ( 11) for several transitions (which, when divided by k B T kin , represents a scaled cross section and has units of Å 2 ). This figure shows the data points computed using the MQCT method in Mandal & Babikov (2023a) and a smooth fitting function from Mandal & Babikov (2023b) that interpolates these data, expands toward the process threshold at low energies, and extrapolates toward the high-energy limit. The resultant values of k n 1 →n ′ 1 (T kin , T rot ) are analyzed below and were used to examine the implications for an astrophysical modeling. T kin and T rot range from 5 to 1000 K.

First, in order to confirm that the new theory and code are correct, we computed the thermal rate coefficients k n 1 →n ′ 1 (T ) at equilibrium, when T kin = T rot , and compared them with the data published previously. In Fig. 2 this comparison is presented for para-H 2 O at three temperatures. In Fig. 3, this information is presented for ortho-H 2 O. The agreement at lower temperatures is reasonably good, but worsens slightly at higher temperatures. The differences for some rate coefficients reach a factor of two to three. This is expected, because in the previously published data, a more approximate method was used to compute the TACS. Namely, in the Eq. ( 13) above, the TACS is a function of two variables, σ n 1 →n ′ 1 (U, T rot ), and the temperature T rot is an input parameter that is independent of U. Importantly, this temperature defines the values of the state populations w n 2 (T rot ) in Eq. ( 13), according to Eq. ( 7). In contrast, in the previous work, we assumed that the average collision energy U(T ) can approximately represent all collisions at a given temperature. This in turn defined the values of the state populations w n 2 (U) and effectively meant that the TACS were a function of one variable only, σ n 1 →n ′ 1 (U), without an explicit dependence on temperature. As a result of these assumptions, in the previous work, the weights w n 2 (U) were different for different energy points U, while here, the weights w n 2 (T rot ) are the same for different energy points because they are detrmined by the chosen rotational temperature, not by the collision energy. This explains some of the differences between the results of this work and previous work shown in Figs. 2 and 3, and it suggests that the approximation employed in previous works is fine at lower temperatures when the distribution of the collision energies is narrow (e.g., at 200 K), but becomes less accurate at higher temperatures when the distribution broadens (e.g., at 800 K). The method implemented here is more accurate, and we recommend following it.

We would also like to note that in our previous work (Mandal & Babikov 2023a,b), we assumed equilibrium between all rotational states of the H 2 O projectile, including the para-and orthostates. This was done largely for historical reasons, following the earlier work on H 2 O + H 2 O by other groups (Buffa et al. 2000;Boursier et al. 2020), which did not distinguish between paraand ortho-states of the H 2 O projectile. Here, we employed a A208, page 4 of 7 rigorous approach by considering para-and ortho-H 2 O projectiles separately, averaging over their states independently, and then computing the weighted average of them using a standard ortho-to-para ratio of 3 (as described in Sect. 2 above). Moreover, we created a user-ready code that allows computing a weighted average of ortho and para projectiles using any other value of the ortho-to-para ratio that is specified as an input parameter. For example, if in a certain model for a certain astrochemical application, we wish to eliminate ortho-H 2 O completely, this can easily be done using the input value of zero for the ortho-to-para ratio. This code is available for download from the GitHubfoot_0 .

In Fig. 4, we present for several of the most intense transitions in H 2 O + H 2 O the dependence of the rate coefficient on both T kin and T rot as a 2D function. In each frame of this figure, the points along the diagonal line correspond to the equilibrium condition, T kin = T rot , while the off-diagonal points describe nonequilibrium situations: T kin < T rot is the area below the diagonal line, and T kin > T rot is the area above the diagonal line. The farther away from the diagonal line, the larger the deviation from equilibrium. The temperatures range from 5 to 1000 K. In many cases, the highest and lowest rate coefficients are observed along the horizontal and vertical axis of the frames in Fig. 4, where one of the two temperatures is relatively low (on the order of 10 K), while the other temperature is far higher (on the order of 100 K A208, page 5 of 7 Mandal, B., et al.: A&A, 688, A208 (2024) line), but for ortho-H 2 O, some of these deviations are larger: 9% for the ground state 1 01 and 22% for the first excited state 1 10 . These differences can be explained by the strong dependence of the excitation rate coefficients on the rotational state of H 2 O.

This comparison highlights the need for accurately determining and constraining the rotational population of the projectiles when modeling a cometary coma. When a thermal distribution of the projectile in the NLTE regime is assumed, this may lead to a significant under-or overestimation of the population of the targets and accordingly, of the abundance of the studied molecules. The case of T rot = 10 K is on the extreme side and could only occur in more distant comets, where T kin is already low, for instance, about 20 K. The other cases considered here (e.g., T rot = 30 K when T kin = 50 K) are more typical, however, and they also demonstrate a non-negligible effect of the NLTE conditions on the population of H 2 O states.

Conclusions

We implemented a two-temperature model to generate a database of rate coefficients for rotational state-to-state transitions in H 2 O + H 2 O collisions that is suitable for the modeling of energy transfer under conditions when the distribution of rotational states of H 2 O deviates from thermodynamic equilibrium. Individual state-to-state transition cross sections from our previous work, predicted using the mixed quantum/classical theory, were re-employed here, but the procedure for averaging over the distribution of rotational states of the projectile H 2 O molecule was modified to permit the deviation of the rotational temperature from the kinetic temperature. Thus, the nascent rate coefficients for rotational transitions in the target H 2 O molecule were obtained as a 2D function of both rotational and translational temperatures.

The results of radiative transfer modeling indicate that for certain conditions that are determined by the column density, significant deviations of the rotational temperature from the kinetic temperature are possible and resultant deviations from the previously used models (which assumed T rot = T kin ) of about 10-20% are typical, but may be even larger for some emission lines of water in the part of the cometary coma where NLTE regime occurs.

These data can be used to model cometary comae and the atmospheres of icy planets where H 2 O + H 2 O collisions are important, but the conditions deviate from thermodynamic equilibrium. For the potential users of this database, we created a computer code that generates rate coefficients for 231 transitions between para-states of water and 210 transitions between orthostates of water (in both excitation and quenching directions) for any given value of kinetic and rotational temperatures in the range from 5 to 1000 K, and for any given value of ortho-to-para ratio of water molecules. This code is available at the GitHubfoot_1 . The data will be submitted to the BASECOL database.