# Rotational state-to-state transition rate coefficients for H₂O + H₂O collisions at nonequilibrium conditions

Bikramaditya Mandal, Michal Zoltowski, Martin Cordiner, François Lique, Dmitri Babikov

► To cite this version:

Bikramaditya Mandal, Michal Zoltowski, Martin Cordiner, François Lique, Dmitri Babikov. Rotational state-to-state transition rate coefficients for H₂O + H₂O collisions at nonequilibrium conditions. Astronomy & Astrophysics - A&A, 2024, 688, 10.1051/0004-6361/202450738. [hal-04699905](https://hal.science/hal-04699905v1)

HAL Id: hal-04699905

https://hal.science/hal-04699905v1

Submitted on 23 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Rotational state-to-state transition rate coefficients for H₂O + H₂O collisions at nonequilibrium conditions

Bikramaditya Mandal<fnref n="1" />, Michal Zoltowski<fnref n="2" />, Martin Cordiner<fnref n="3" />,<fnref n="4" />, Francois Lique<fnref n="5" />, and Dmitri Babikov<fnref n="1" />

Received 16 May 2024 / Accepted 5 July 2024

## ABSTRACT

**Aims.** The goal is to develop a database of rate coefficients for rotational state-to-state transitions in H₂O + H₂O collisions that is suitable for the modeling of energy transfer in nonequilibrium conditions, in which the distribution of rotational states of H₂O deviates from local thermodynamic equilibrium.

**Methods.** A two-temperature model was employed that assumed 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<sub>kin</sub> and T<sub>rot</sub>.

**Results.** Upon introducing our new parameterization of the collisional rates, taking into account their dependence on both T<sub>kin</sub> and T<sub>rot</sub>, we find a change of up to 20% in the H₂O rotational level populations for both ortho and para-H₂O for the part of the cometary coma where the nonequilibrium regime occurs.

**Key words.** molecular data – molecular processes – radiative transfer – comets: general

## 1. Introduction

Inelastic collisions between water molecules, H₂O + H₂O, play an important role in collisional energy transfer in environments such as atmospheres of icy planets (Wiström et al. 2020; Vorburger et al. 2022), in cometary comae (Cochran et al. 2015; Cordiner et al. 2022), and in the atmospheres of water-rich 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₂O + H₂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₂O + H₂O system as a whole. Recently, noticeable progress in the modeling of H₂O + H₂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₂O and 210 transitions in ortho-H₂O due to collisions with a thermalized bath of H₂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₂O + H₂O system involve averaging over the rotational states of the projectile H₂O assuming a Boltzman distribution of these states at given kinetic temperature, that is, the equilibrium conditions. These data can probably

<footnote n="1">¹ Chemistry Department, Wehr Chemistry Building, Marquette University, Milwaukee, Wisconsin 53201-1881, USA
e-mail: dmitri.babikov@marquette.edu</footnote>

<footnote n="2">² Institute of Physics, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziącka Street 5,
87-100 Toruń, Poland</footnote>

<footnote n="3">³ Solar System Exploration Division, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="4">⁴ Institute for Astrophysics and Computational Sciences, The Catholic University of America, Washington, DC 20064, USA</footnote>

<footnote n="5">⁵ University of Rennes, CNRS, IPR (Institut de Physique de Rennes) - UMR 6251, 35000 Rennes, France</footnote>
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₂O + H₂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₂O + H₂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₂O). Then, in the modeling of radiation transfer, these rate coefficients might be averaged over the states of the projectile H₂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₂O + H₂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₂O molecules, while $T_{rot}$ can be used to average over the rotational states of the projectile H₂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₂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₃O⁺, H₂O, and HDO collided with H₂ 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₂O + H₂O, and we explore the implications of nonequilibrium conditions for an astrophysical modeling of H₂O-rich environments, such as the atmospheres of icy planets and cometary comae.

## 2. Details of the method

The standard procedure for predicting thermal state-to-state transition rate coefficients $k_{n_1 \to n'_1}(T)$ begins with averaging the state-to-state transition cross sections $\sigma_{n_1 n_2 \to n'_1 n'_2}(E)$ over the Boltzmann-Maxwell distribution of collision energies $E$ at a certain temperature $T$,

$$k_{n_1 n_2 \to n'_1 n'_2}(T) = \frac{v_{\text{ave}}(T)}{(k_B T)^2} \int_{E=0}^{\infty} E \sigma_{n_1 n_2 \to n'_1 n'_2}(E) e^{-\frac{E}{k_B T}} dE. \quad (1)$$

In this expression, $v_{\text{ave}}(T) = \sqrt{8k_B T/\pi\mu}$ is the average collision velocity, $k_B$ is the Boltzmann constant, $\mu$ 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-top-rotor 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,

$$k_{n_1 \to n'_1}(T) = \sum_{n_2} w_{n_2}(T) \sum_{n'_2} k_{n_1 n_2 \to n'_1 n'_2}(T), \quad (2)$$

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

$$w_{n_2}(T) = \frac{(2j_2 + 1)e^{-\frac{E_2}{k_B T}}}{Q_2(T)}. \quad (3)$$

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,

$$Q_2(T) = \sum_{n_2} (2j_2 + 1) e^{-\frac{E_2}{k_B T}}. \quad (4)$$

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 ortho-states of the projectile as if they would be different species, which gives two values of $k_{n_1 \to n'_1}(T)$ for each $n_1 \to 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₂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 \to 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,

$$k_{n_1 n_2 \to n'_1 n'_2}(T_{\text{kin}}) = \frac{v_{\text{ave}}(T_{\text{kin}})}{(k_B T_{\text{kin}})^2} \int_{E=0}^{\infty} E \sigma_{n_1 n_2 \to n'_1 n'_2}(E) e^{-\frac{E}{k_B T_{\text{kin}}}} dE. \quad (5)$$

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

$$Q_2(T_{\text{rot}}) = \sum_{n_2} (2j_2 + 1) e^{-\frac{E_2}{k_B T_{\text{rot}}}}, \quad (6)$$

$$w_{n_2}(T_{\text{rot}}) = \frac{(2j_2 + 1) e^{-\frac{E_2}{k_B T_{\text{rot}}}}}{Q_2(T_{\text{rot}})}, \quad (7)$$

and is used for averaging over this distribution,

$$k_{n_1 \to n'_1}(T_{\text{kin}}, T_{\text{rot}}) = \sum_{n_2} w_{n_2}(T_{\text{rot}}) \sum_{n'_2} k_{n_1 n_2 \to n'_1 n'_2}(T_{\text{kin}}). \quad (8)$$

In this way, the value of state-to-state transition rate coefficient depends on two temperatures, $T_{kin}$ and $T_{rot}$. At equilibrium, when
$T_{\text{rot}} = T_{\text{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₂O + H₂O. The reason is the huge number of individual state-to-state transitions $n_1 n_2 \rightarrow 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 ortho-ortho lower-energy transitions in the target H₂O, and we used a rotational basis set of 38 para- and 38 ortho-states for the projectile H₂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 $\sigma_{n_1 n_2 \rightarrow 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

$$k_{n_1 \rightarrow n'_1}(T_{\text{kin}}, T_{\text{rot}}) = \frac{v_{\text{ave}}(T_{\text{kin}})}{(k_B T_{\text{kin}})^2} \int_{E=0}^{\infty} \sigma_{n_1 \rightarrow n'_1}(E, T_{\text{rot}}) e^{-\frac{E}{k_B T_{\text{kin}}}} dE. \quad (9)$$

We introduced a thermally averaged cross section (TACS) for the transition $n_1 \rightarrow n'_1$, defined as

$$\sigma_{n_1 \rightarrow n'_1}(E, T_{\text{rot}}) = \sum_{n_2} w_{n_2}(T_{\text{rot}}) \sum_{n'_2} \sigma_{n_1 n_2 \rightarrow n'_1 n'_2}(E). \quad (10)$$

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_{\text{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

<figure><img src="image_2.png" /><figcaption>Fig. 1. Function F, which includes all factors under the integral of Eq. (11) for several transitions between the ground state $0_{00}$ of para-H₂O and its several low-lying rotational states (as indicated by different colors in the figure) at $T_{\text{kin}} = 100$ K. The solid, dashed, and dash-dotted lines correspond to a rotational temperature of $T_{\text{rot}} = 10, 50, \text{ and } 100$ K, respectively.</figcaption></figure>

(Mandal & Babikov 2023b),

$$k_{n_1 \rightarrow n'_1}(T_{\text{kin}}, T_{\text{rot}}) = \frac{v_{\text{ave}}(T_{\text{kin}})}{(k_B T_{\text{kin}})^2} \left( \frac{e^{-\frac{\Delta E}{2k_B T_{\text{kin}}}}}{2j + 1} \right) \times \int_{U_{\text{min}}}^{\infty} \tilde{\sigma}_{n_1 n'_1}(U, T_{\text{rot}}) e^{-\frac{U}{k_B T_{\text{kin}}}} [1 + (\frac{\Delta E}{4U})^2] [1 - (\frac{\Delta E}{4U})^2] dU. \quad (11)$$

Here, U is the kinetic energy of the mixed quantum/classical trajectories, and $U_{\text{min}} = \Delta E/4$ is its threshold value for a given $n_1 \rightarrow n'_1$ transition. $\Delta 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),

$$\tilde{\sigma}_{n_1 n'_1}(U, T_{\text{rot}}) = \frac{1}{2}[(2j + 1)\sigma_{n_1 \rightarrow n'_1}(U, T_{\text{rot}}) + (2j' + 1)\sigma_{n'_1 \rightarrow n_1}(U, T_{\text{rot}})]. \quad (12)$$

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,

$$\sigma_{n_1 \rightarrow n'_1}(U, T_{\text{rot}}) = \sum_{n_2} w_{n_2}(T_{\text{rot}}) \sum_{n'_2} \sigma_{n_1 n_2 \rightarrow n'_1 n'_2}(U). \quad (13)$$

## 3. Results and discussion

From previous work (Mandal & Babikov 2023a), the cross-sections $\sigma_{n_1 n_2 \rightarrow n'_1 n'_2}$ are available at six collision energies: $U = 133, 200, 267, 400, 533, \text{ and } 708 \text{ cm}^{-1}$. These were used to compute six values of $\sigma_{n_1 \rightarrow n'_1}(U, T_{\text{rot}})$ at any given input temperature
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Comparison of the rate coefficients computed here for the case of thermodynamic equilibrium $T_{\text{kin}} = T_{\text{rot}}$ (vertical) vs. thermal rate coefficients reported previously (Mandal & Babikov 2023b; horizontal) for the quenching of 231 states of para-H₂O at three temperatures.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Same as Fig. 2, but for the quenching of 210 states of ortho-H₂O.</figcaption></figure>

$T_{\text{rot}}$ using Eq. (10). The excitation and quenching cross sections were combined in Eq. (12) to give $\bar{\sigma}_{n_1 n'_1}(U, T_{\text{rot}})$. These were interpolated and integrated over $U$ for any given input temperature $T_{\text{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_{\text{kin}}$, represents a scaled cross section and has units of $\AA^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 \to n'_1}(T_{\text{kin}}, T_{\text{rot}})$ are analyzed below and were used to examine the implications for an astrophysical modeling. $T_{\text{kin}}$ and $T_{\text{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 \to n'_1}(T)$ at equilibrium, when $T_{\text{kin}} = T_{\text{rot}}$, and compared them with the data published previously. In Fig. 2 this comparison is presented for para-H₂O at three temperatures. In Fig. 3, this information is presented for ortho-H₂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, $\sigma_{n_1 \to n'_1}(U, T_{\text{rot}})$, and the temperature $T_{\text{rot}}$ is an input

parameter that is independent of $U$. Importantly, this temperature defines the values of the state populations $w_{n_2}(T_{\text{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, $\sigma_{n_1 \to 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_{\text{rot}})$ are the same for different energy points because they are determined 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₂O projectile, including the para- and ortho-states. This was done largely for historical reasons, following the earlier work on H₂O + H₂O by other groups (Buffa et al. 2000; Boursier et al. 2020), which did not distinguish between para- and ortho-states of the H₂O projectile. Here, we employed a
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** 2D dependence of the nonequilibrium rate coefficients for 16 strong transitions in the H₂O + H₂O system, as indicated in the figure. The rotational and kinetic temperatures are plotted along the horizontal and vertical axis, respectively. Color is used to show the magnitude of rate coefficients in units of cm³ s⁻<fnref n="1" /> in log scale as indicated in the color bar.</figcaption></figure>

rigorous approach by considering para- and ortho-H₂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₂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 GitHub<fnref n="1" />.

In Fig. 4, we present for several of the most intense transitions in H₂O + H₂O the dependence of the rate coefficient on both T<sub>kin</sub> and T<sub>rot</sub> as a 2D function. In each frame of this figure, the points along the diagonal line correspond to the equilibrium condition, T<sub>kin</sub> = T<sub>rot</sub>, while the off-diagonal points describe nonequilibrium situations: T<sub>kin</sub> < T<sub>rot</sub> is the area below the diagonal line, and T<sub>kin</sub> > T<sub>rot</sub> 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).

<footnote n="1">¹ https://github.com/MarquetteQuantum/MOLRATES</footnote>
and above). We found that in extreme nonequilibrium conditions like this, the rate coefficients for many transitions can deviate significantly from their equilibrium values.

## 4. Nonequilibrium radiative transfer modeling

In order to determine the impact of the rotational temperature of H₂O projectiles on the excitation conditions of H₂O targets in nonequilibrium conditions, we performed radiative transfer modeling using the new collisional data.

The coupled radiative transfer and statistical equilibrium equations were solved using the escape probability approximation at steady state, implemented in the public version of the RADEX code (Van der Tak et al. 2007). We selected physical conditions corresponding to a typical comet located at $R_h = 1$ AU from the Sun with a total production rate from the nucleus $Q_{H2O} = 10^{29}$ s⁻¹, an expansion velocity $v_{exp} = 0.8$ km s⁻¹, and at two different gas kinetic temperatures, $T_{kin} = 50$ and 100 K. Water was the only projectile considered, and the H₂O density was described by a standard isotropic Haser model (Haser 1957). The H₂O column density was set to $10^{17}$ cm⁻², which corresponds to the typical column density in the part of the coma where NLTE effects are expected to be important. At this column density, the opacities of the H₂O lines will be high. We then anticipate that the transition from LTE to NLTE will occur at a lower H₂O density than in the optically thin regime. For the background radiation, we incorporated in the calculations the cosmic microwave radiation (2.73 K) along with solar radiation. Regarding the latter, we adopted the same approach as Faure et al. (2020) and used a diluted blackbody radiation of 5770 K with a dilution factor $W = \Omega_s/4\pi$, where $\Omega_s$ equals $6.79 \times 10^{-5}$ sr, representing the solid angle of the Sun as seen from 1 AU.

The results of our NLTE models are presented in Figs. 5 and 6, where the populations of several lower-energy rotational levels of H₂O (both para- and ortho-H₂O) are plotted as a function of the H₂O density for $T_{kin} = 50$ and 100 K. We selected these levels since they are the most populated levels in typical physical coma conditions, but also because they are all radiatively connected to one or more rotational states and then would be involved in radiative transitions detected by telescopes. The distribution of the H₂O rotational levels of the projectiles was performed assuming a rotational temperature of $T_{rot} = 10, 30, \text{ and } 50$ K at $T_{kin} = 50$ K and of $T_{rot} = 10, 50, \text{ and } 100$ K at $T_{kin} = 100$ K. When $T_{rot}$ and $T_{kin}$ are equal, this is equivalent to the thermal equilibrium distribution of the rotational states of the projectiles.

First, three different regimes exist: the fluorescence equilibrium at low densities ($n_{H2O} < 10$ cm⁻³), the NLTE regime in the range $10 < n_{H2O} < 10^6$ cm⁻³, and the LTE regime at higher densities. Therefore, the NLTE regime is expected for the intermediate part of the coma, and this conclusion is relatively independent of the rotational distribution of the projectiles. We also note that the energy level populations at fluorescence equilibrium do not significantly vary with temperature. In contrast, the LTE population distribution is sensitive to the temperature.

Then, in the intermediate (NLTE) regime, the population of the energy levels is dependent on the rotational temperature of the projectile. At a lower kinetic temperature ($T_{kin} = 50$ K, Fig. 5), we found that for para-H₂O, large deviations exist for the energy level populations in this regime (rotational temperature of 10 K vs. thermalized projectiles). Namely, near $n_{H2O} = 10^3$ cm⁻³, the deviation of the state populations (solid vs. dash-dotted line) reaches 21% for the ground state $\Omega_{00}$ and 54% for the

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Relative populations of several low-lying rotational states of para- and ortho-H₂O (targets) as a function of the H₂O (projectile) density at $T_{kin} = 50$ K. The solid, dashed, and dash-dotted lines correspond to a rotational temperature of $T_{rot} = 10, 30, \text{ and } 50$ K, respectively.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Same as in Fig. 5, but at $T_{kin} = 100$ K and $T_{rot} = 10, 50, \text{ and } 100$ K represented by solid, dashed, and dash-dotted lines, respectively.</figcaption></figure>

first excited state $|1_{11}\rangle$. In contrast, for ortho-H₂O, these deviations are smaller, about 3 and 9% for states $|0_{11}$ and $|1_{10}$, respectively. However, as the temperature increases, their behavior reverses. At $T_{kin} = 100$ K (in Fig. 6), we found that when the rotational temperature is 10 K, the state populations for para-H₂O are within 10% of the thermalized case (again, solid vs. dash-dotted
line), but for ortho-H₂O, some of these deviations are larger: 9% for the ground state l₀₁ and 22% for the first excited state l₁₀. These differences can be explained by the strong dependence of the excitation rate coefficients on the rotational state of H₂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<sub>rot</sub> = 10 K is on the extreme side and could only occur in more distant comets, where T<sub>kin</sub> is already low, for instance, about 20 K. The other cases considered here (e.g., T<sub>rot</sub> = 30 K when T<sub>kin</sub> = 50 K) are more typical, however, and they also demonstrate a non-negligible effect of the NLTE conditions on the population of H₂O states.

**5. Conclusions**

We implemented a two-temperature model to generate a database of rate coefficients for rotational state-to-state transitions in H₂O + H₂O collisions that is suitable for the modeling of energy transfer under conditions when the distribution of rotational states of H₂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₂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₂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<sub>rot</sub> = T<sub>kin</sub>) 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₂O + H₂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 ortho-states 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 GitHub<fnref n="2" />. The data will be submitted to the BASECOL database.

*Acknowledgements.* This research was supported by NASA Astrophysics Program, grant number 80NSSC24K0208, and by French INSU-CNRS "programme national de planétologie, PNP". D.B. acknowledges the support of Way Klingler Fellowship and Haberman research fund. F.L. acknowledges Rennes Metropole and the Institut Universitaire de France for financial support. M.A.C. was supported by the National Science Foundation under Grant No. AST-2009253. We used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-5CH11231. This research also used HPC resources at Marquette funded in part by the National Science foundation award CNS-1828649.

**References**

Agg, P., & Clary, D. 1991, J. Chem. Phys., 95, 1037

Alexander, M., Dagdigian, P. J., Werner, H.-J., et al. 2023, Comp. Phys. Commun., 289, 108761

Boursier, C., Mandal, B., Babikov, D., & Dubernet, M. 2020, MNRAS, 498, 5489

Brinch, C., & Hogerheijde, M. 2010, A&A, 523, A25

Buffa, G., Tarrini, O., Scappini, F., & Cecchi-Pestellini, C. 2000, ApJS, 128, 597

Cochran, A. L., Levasseur-Regourd, A.-C., Cordiner, M., et al. 2015, Space Sci. Rev., 197, 9

Cordiner, M., Coulson, I., Garcia-Berrios, E., et al. 2022, ApJ, 929, 38

Demes, S., Lique, F., Faure, A., & van der Tak, F. F. 2023, MNRAS, 518, 3593

Dubernet, M., Boursier, C., Denis-Alpizar, O., et al. 2024, A&A, 683, A40

Faure, A., Lique, F., & Loreau, J. 2020, MNRAS, 493, 776

Faure, A., Žóltowski, M., Wiesenfeld, L., Lique, F., & Bergeat, A. 2024, MNRAS, 527, 3087

García Muñoz, A., Ramos, A. A., & Faure, A. 2024, Icarus, 415, 116080

Haser, L. 1957, Bulletins de l'Académie Royale de Belgique, 43, 740

Hutson, J. M., & Le Sueur, C. R. 2019, Comput. Phys. Commun., 241, 9

Hutson, J. M., & Le Sueur, C. R. 2020, MOLSCAT, Bound and Field, version 2020.0, https://github.com/molscat/molscat

Lique, F., Honvault, P., & Faure, A. 2014, Int. Rev. Phys. Chem., 33, 125

Mandal, B., & Babikov, D. 2023a, A&A, 671, A51

Mandal, B., & Babikov, D. 2023b, A&A, 678, A51

Mandal, B., Joy, C., Bostan, D., Eng, A., & Babikov, D. 2023, J. Phys. Chem. Lett., 14, 817

Mandal, B., Bostan, D., Joy, C., & Babikov, D. 2024, Comput. Phys. Commun., 294, 108938

Raizer, Y. P., & Allen, J. E. 1997, Gas Discharge Physics (Berlin: Springer), 2

Ramos, A. A., & Elitzur, M. 2018, A&A, 616, A131

Van der Tak, F., Black, J. H., Schöier, F., Jansen, D., & van Dishoeck, E. F. 2007, A&A, 468, 627

van der Tak, F. F., Lique, F., Faure, A., Black, J. H., & van Dishoeck, E. F. 2020, Atoms, 8, 15

Vorburger, A., Fatemi, S., Galli, A., et al. 2022, Icarus, 375, 114810

Wirstöm, E., Bjerkeli, P., Rezac, L., Brinch, C., & Hartogh, P. 2020, A&A, 637, A90

<footnote n="2">² https://github.com/MarquetteQuantum/MOLRATES</footnote>