NLTE modelling of water-rich exoplanet atmospheres.
Cooling and heating rates

A. García Muñoz, A. Asensio Ramos, Alexandre Faure

► To cite this version:

A. García Muñoz, A. Asensio Ramos, Alexandre Faure. NLTE modelling of water-rich exoplanet atmospheres. Cooling and heating rates. Icarus, 2024, 415, pp.116080. 10.1016/j.icarus.2024.116080. hal-04738719

HAL Id: hal-04738719

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

Submitted on 17 Jul 2025

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 Attribution 4.0 International License
NLTE modelling of water-rich exoplanet atmospheres. Cooling and heating rates.

A. GARCÍA MUÑOZ,¹ A. ASENSIO RAMOS,² AND A. FAURE³

¹ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France

² Instituto de Astrofísica de Canarias, 38205, La Laguna, Tenerife, Spain; Departamento de Astrofísica, Universidad de La Laguna, 38205, La Laguna, Tenerife, Spain

³ Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

(Accepted April 4, 2024)

Submitted to Icarus

ABSTRACT

The hydrogen and water molecules respond very differently to the collisional-radiative processes taking place in planetary atmospheres. Naturally, the question arises whether H₂O-rich atmospheres are more (or less) resilient to long-term mass loss than H₂-dominated ones if they radiate away the incident stellar energy more (or less) efficiently. If confirmed, the finding would have implications on our understanding of the evolution of sub-Neptune exoplanets. As a key step towards answering this question, we present a non-local thermodynamic equilibrium (NLTE) model of H₂O for the atmospheric region where the gas accelerates to escape the planet and conditions relevant to close-in sub-Neptunes. Our exploratory calculations for isothermal gas composed of H₂, H₂O and e⁻ reveal that: 1) In the pressure region ~10⁻²–10⁻⁴ dyn cm⁻² where the stellar extreme-ultraviolet (XUV) photons are typically deposited in the atmosphere, H₂O is in rotational LTE but vibrational NLTE. Vibrational LTE is facilitated by high H₂O abundances and fractional ionizations, and we report critical densities for the LTE-NLTE transition; 2) Vibrational cooling may locally dominate over rotational cooling, partly because of the comparatively small opacities of ro-vibrational lines; 3) Even low H₂O abundances notably enhance the cooling, foreseeably offsetting some of the stellar heating; 4) Heating due to the deposition of stellar infrared (IR) photons is significant at pressures ≳0.1 dyn cm⁻². We estimate the contribution of H₂O excitation to the internal energy of the gas and speculate on the photodissociation from the excited vibrational states. Ultimately, our findings motivate the consideration of NLTE in the mass loss rate calculations of H₂O-rich atmospheres.

Keywords: ...

Corresponding author: Antonio García Muñoz
antonio.garciamunoz@cea.fr, tonhingm@gmail.com
# 1. INTRODUCTION

Water is thought to be abundant at some low-mass exoplanets (Venturini et al. 2020; Bitsch et al. 2021; Kite & Schaefer 2021; Kimura & Ikoma 2022; Izidoro & Piani 2023). Proving the reality of H₂O-rich planets is challenging though as mass and size alone, two properties accurately known for many planets, are often insufficient to conclude on their interior structures and atmospheric compositions (Valencia et al. 2007; Acuña et al. 2021; Delrez et al. 2021). Even when the transmission spectrum shows evidence for H₂O, the intrinsic degeneracies in the interpretation of the spectrum make it difficult to decide whether H₂O is a major or minor atmospheric constituent (Benneke et al. 2019; Madhusudhan et al. 2020). The long-term stability of the atmosphere provides a complementary way of assessing its composition. There is indeed a growing sample of planets for which the existence of atmospheres is firmly established but they are unlikely to be H₂-dominated as they would have been lost to space long ago under their current irradiation conditions. That leaves H₂O and other gases of high mean molecular weight as more plausible alternatives (García Muñoz et al. 2020, 2021; Piaulet et al. 2023).

The community is in the process of testing these and other ideas about the atmospheric composition of low-mass planets with the James Webb Space Telescope (JWST) (Ih et al. 2023; Lincowski et al. 2023; Moran et al. 2023) and other facilities. Simultaneously, we must assess the validity of the assumptions made in the models used for predicting the mass loss rates and how they affect the stability of the atmospheres. For example, a H₂O-rich atmosphere is expected to respond to infrared (IR) radiation very differently than a H₂-dominated one because, unlike the H₂ molecule, H₂O has a large dipole moment. Also, H₂O has numerous low-energy states that are easily populated in collisions with molecules (including atoms) and electrons and that, once excited, radiate promptly. In particular, self-collisions with other H₂O molecules are efficient at populating the H₂O states.

Our work assesses the significance of bound-bound collisional-radiative processes within the H₂O molecule for the cooling and heating (net cooling, hereafter) of H₂O-rich atmospheres. We build for that purpose a model that solves the Non-Local Thermodynamic Equilibrium (NLTE) problem of H₂O under quite general conditions. The NLTE treatment is warranted at the low pressures where the atmospheric gas is accelerated outwards to escape the planet and the population of H₂O states may more easily depart from a Boltzmann distribution. The model considers the collisions of H₂O with the background gas, and the interaction with upwelling radiation emitted from deeper atmospheric layers, with internally-generated diffuse radiation and with stellar radiation. The radiative interactions occur in the far-IR for purely rotational transitions (that occur with no change in the vibrational quantum numbers) and in the near- and mid-IR for ro-vibrational transitions (that occur with changes in the vibrational quantum numbers). To our knowledge, a similar investigation of exoplanet atmospheres has never been done. We are exploring the implications on various atmospheric properties and on the demographics of sub-Neptunes in ongoing work. The connected nature of the energy budget, chemistry and dynamics of atmospheres (Sánchez-Lavega et al. 2023) suggests that H₂O NLTE might alter the three-dimensional structure of exoplanet atmospheres, especially at low and intermediate pressures. The idea remains to be investigated by existent General Circulation Models (Kaspi & Showman 2015; Wordsworth 2015; Carone et al. 2018).

The H₂O NLTE problem is certainly not new in the planetary sciences and astrophysics. It has been treated in the remote sensing of the Earth's atmosphere to constrain its temperature and energy budget and to infer H₂O abundances (López-Puertas & Taylor 2001; Feofilov et al. 2009). It has been considered in interstellar and circumstellar clouds, stars and galaxies (Cernicharo & Crovisier 2005; van Dishoeck et al. 2013; Gray et al. 2016), prompting various calculations of cooling rates (Goldsmith & Langer 1978; Hollenbach & McKee 1979; Neufeld & Kaufman 1993; Morris et al. 2009). In the latter astrophysical environments, H₂O occurs in trace amounts and its collisional excitation is dominated by H₂ and, possibly, H atoms or electrons. The H₂O NLTE problem has also been treated in the investigation of the envelopes of comets and icy moons (Crovisier 1984; Xie & Mumma 1992; Gersch et al. 2018; Villanueva et al. 2023), where H₂O is a major constituent. These works often focus on predicting the populations of excited states for comparison against remote sensing measurements. The problem addressed here is somewhat different and more general. We consider a broad range of compositions, from H₂O being a trace constituent to being dominant, and the possibility of both the purely rotational and ro-vibrational excitation of H₂O in collisions with H₂, H₂O and electrons of kinetic temperatures from 200 to 1,500 K. Our work focuses on the cooling rates rather than on the specifics of the H₂O spectrum.
## 2. THE H₂O MOLECULE

H₂O occurs in ortho (o-H₂O; total nuclear spin of the hydrogen atoms $I=1$; nuclear spin statistical weight $g_I=3$) and para (p-H₂O; $I=0$; $g_I=1$) isomers. The two isomers do not readily interconvert through collisions or radiation, and we treat them separately. As a non-linear triatomic molecule, H₂O possesses 3 vibrational modes of motion, identified by the quantum numbers $v=(v₁v₂v₃)$ for symmetric stretching ($v₁$), bending ($v₂$), and asymmetric stretching ($v₃$). The rotational motion is described by the quantum numbers ($J K_a K_c$), or alternatively ($J_τ$ with $τ=K_a-K_c$). $J$ is the total rotational quantum number, and $K_a$ and $K_c$ are the respective projections of $J$ on the molecular axes with the smallest and greatest moments of inertia.

Building the NLTE model requires specifying the participating o- and p-H₂O states, the collisional-radiative processes through which they interact with the medium and their rate coefficients. Full details on our molecular model are given in Appendix A. In short, it consists for o-H₂O of ~400 states within the five lowest-energy vibrational states, namely $v=(000)$, (010), (020), (100) and (001). They are connected through ~7,600 radiative transitions. It considers the collisions of o-H₂O with the molecules H₂ and H₂O and with electrons. About 83,000 channels for inelastic collisions of o-H₂O with electrons are included, and comparable numbers for collisions with H₂ and H₂O. Inelastic collisions with ions are omitted but this should not be a limitation of the model as the ion densities in the atmosphere are small where the H₂O molecule remains undissociated, and inelastic collisions with electrons are dominant once ionization becomes significant. The number of states, radiative transitions and collisional channels for p-H₂O are similar to those for o-H₂O.

## 3. COOLING RATE

We review the formulation of the NLTE problem and its contribution to the net cooling of the gas. The treatment is standard and we include it here for later reference.

We are concerned with the energy equation for a radiating fluid (Zel'dovich & Raizer 2002):

$$ \frac{\partial \rho E}{\partial t} + \nabla \cdot [\rho H \mathbf{u} + \mathbf{q}] - \rho \mathbf{f}_{\text{ext}} \cdot \mathbf{u} = -\Gamma = - \int \frac{dI_{\lambda}}{ds} d\Omega d\lambda. \quad (1) $$

Here, $\rho E$ is the total energy, including the bulk kinetic energy and the internal energy; $\rho H=\rho E + p$ is the total enthalpy, and $p$ is the pressure; $\mathbf{u}$ is the bulk velocity of the gas, and $\rho \mathbf{f}_{\text{ext}} \cdot \mathbf{u}$ is the work of the external forces; $\mathbf{q}$ is the heat flux. The internal energy contains the enthalpy of formation (plus dissociation and ionization) of the different particles (molecules and electrons), and their translational, rotational, vibrational and electronic energies (as applicable).

The right hand side of Eq. 1 represents the exchange of energy between the gas and its surroundings through radiation. $\Gamma$ is calculated by solving the radiative transfer equation:

$$ \frac{dI_{\lambda}}{ds} = -\kappa_{\lambda} I_{\lambda} + \varepsilon_{\lambda}, \quad (2) $$

which expresses the variation in radiance $I_\lambda$ per unit length $ds$, for absorption and emission coefficients $\kappa_\lambda$ and $\varepsilon_\lambda$, and integrating it over solid angle $\Omega$ and wavelength $\lambda$. Radiation occurs over bound-bound, bound-free/free-bound and free-free transitions, and $\kappa_\lambda=\kappa_\lambda^{\text{BB}}+\kappa_\lambda^{\text{BF}}+\kappa_\lambda^{\text{FF}}$ and $\varepsilon_\lambda=\varepsilon_\lambda^{\text{BB}}+\varepsilon_\lambda^{\text{BF}}+\varepsilon_\lambda^{\text{FF}}$, such that:

$$ \frac{dI_{\lambda}}{ds} = \underbrace{-\kappa_{\lambda}^{\text{BB}} I_{\lambda} + \varepsilon_{\lambda}^{\text{BB}}}_{\rightarrow \Gamma^{\text{BB}}} - \underbrace{-\kappa_{\lambda}^{\text{BF}} I_{\lambda} + \varepsilon_{\lambda}^{\text{BF}}}_{\rightarrow \Gamma^{\text{BF/FB}}} - \underbrace{-\kappa_{\lambda}^{\text{FF}} I_{\lambda} + \varepsilon_{\lambda}^{\text{FF}}}_{\rightarrow \Gamma^{\text{FF}}}, \quad (3) $$

and separate contributions $\Gamma^{\text{BB}}$, $\Gamma^{\text{BF/FB}}$ and $\Gamma^{\text{FF}}$ can be defined. In practice, the solution to the radiative transfer problem for the bound-bound transitions can be decoupled from the general radiative transfer problem if at the relevant wavelengths $\kappa_\lambda \approx \kappa_\lambda^{\text{BB}}$ and $\varepsilon_\lambda \approx \varepsilon_\lambda^{\text{BB}}$. This is the case for the problem of interest here because the bound-bound radiative transitions of H₂O occur mostly in the IR, whereas the bound-free/free-bound transitions occur in the ultraviolet and the free-free transitions typically contribute weakly across the broad spectrum. $\Gamma^{\text{BF/FB}}$ and $\Gamma^{\text{FF}}$ are
often incorporated in the hydrodynamical models that solve Eq. 1, whereas $\Gamma^{BB}$ is omitted or severely simplified in them. Our focus in this work is to calculate $\Gamma^{BB}$ for a variety of conditions relevant to H$_2$O-rich atmospheres.

The absorption and emission coefficients for bound-bound transitions are given by sums over all the H$_2$O states of indices i, j, (number) densities $n_i$, $n_j$ and energies $E_j>E_i$ that are connected radiatively:

$$ \kappa_{\lambda}^{BB} = \sum_{i,j>i} \frac{h\lambda_{ij}}{4\pi} (n_i B_{ij} - n_j B_{ji}) \psi_{\lambda}^{ij} \quad (4) $$

$$ \varepsilon_{\lambda}^{BB} = \sum_{i,j>i} \frac{hc}{4\pi\lambda_{ij}} n_j A_{ji} \psi_{\lambda}^{ij} \quad (5) $$

where h and c are the Planck constant and speed of light, respectively, and $\lambda_{ij}=hc/(E_j-E_i)$ is the wavelength. $A_{ji}$, $B_{ji}$ and $B_{ij}$ are the Einstein coefficients for spontaneous emission, induced emission and absorption, related through $B_{ji}/A_{ji}=\lambda_{ij}^3/2hc$ and $g_jB_{ji}=g_iB_{ij}$. The above expressions assume complete frequency redistribution, with absorption and emission described by the same line profile $\psi_{\lambda}^{ij}$.

Following the above, the net cooling rate for bound-bound transitions, or net cooling rate for short, is given by:

$$ \Gamma^{BB} = \sum_{i,j>i} \left[ n_j A_{ji} - (n_i B_{ij} - n_j B_{ji}) \Phi_{ij} \frac{\lambda_{ij}^2}{c} \right] (E_j - E_i), \quad (6) $$

where:

$$ \Phi_{ij} = \frac{1}{4\pi} \iint I_{\lambda} \psi_{\lambda}^{ij} d\Omega d\lambda $$

is the mean line intensity integrated over wavelength. Alternatively:

$$ \Gamma^{BB} = \sum_{i,j>i} n_j A_{ji} p_{ji} (E_j - E_i), \quad (7) $$

and the so-called net radiative bracket (Athay & Skumanich 1971; Elitzur & Asensio Ramos 2006):

$$ p_{ji} = 1 - \frac{n_i g_j - n_j g_i \lambda_{ij}^5 \Phi_{ij}}{n_j g_i 2hc^2}. \quad (8) $$

As used here, $p_{ji}$ includes in $\Phi_{ij}$ the contributions from both external radiation entering the atmosphere and diffuse radiation generated within. We note however that the usual definition of $p_{ji}$ includes only the contribution from diffuse radiation (Elitzur & Asensio Ramos 2006), in which case it is related to the probability for the line photon to escape the gas (Athay & Skumanich 1971; Irons 1978), whereas the external contribution is treated through a separate term. This being a matter of nomenclature that affects only the calculations for which the gas is externally irradiated, we have followed the above definition because the corresponding $p_{ji}$ are easy to extract from our NLTE solver and because keeping separate terms does not add any physical insight to the discussion. It is sometimes useful to separate $\Gamma^{BB}$ into the components $\Gamma_v^{BB}$, each of them collecting the cooling rate from radiative transitions with upper states in the specified vibrational state v. By construction, the summation of the rotational component $\Gamma_{000}^{BB}$ and the ro-vibrational components $\Gamma_{010}^{BB}$, $\Gamma_{020}^{BB}$, $\Gamma_{100}^{BB}$ and $\Gamma_{001}^{BB}$ amounts to the net cooling rate.

Figure 1 sketches the H$_2$O molecular model. The system is assumed closed, meaning that any collisional-radiative process that initiates from a bound state i results into another bound state j. This generally requires that the chemical processes producing or destroying the H$_2$O molecule proceed more slowly than the relaxation within it, thereby rendering the two problems separable. This condition is approximately met in the region of an atmosphere where H$_2$O remains undissociated and where the relaxation times are shorter than the chemical loss times. For our NLTE calculations, we have assumed that the total density of H$_2$O is known and satisfies $[H_2O]=\sum_i n_i$.
The density of a state i is determined by solving the steady-state problem:

$$ \frac{dn_i}{dt} = 0 = \sum_{j>i} n_j A_{ji} p_{ji} - \sum_{j<i} n_i A_{ij} p_{ij} + \sum_{j \neq i} (n_j C_{ji} - n_i C_{ij}). \quad (9) $$

for the balance of production and loss rates through all possible collisional-radiative processes. $C_{ji}$ is the frequency for collisional deexcitation from j to i for $j>i$, and the frequency for collisional excitation from j to i for $j<i$. Summation of this equation $\times E_i$ over all the states results in the equivalent form for the net cooling rate:

$$ \Gamma^{\text{BB}} = \sum_{i,j>i} (n_i C_{ij} - n_j C_{ji}) (E_j - E_i). \quad (10) $$

Cooling ($\Gamma^{\text{BB}} > 0$) occurs when excitation dominates over deexcitation, and heating ($\Gamma^{\text{BB}} < 0$) when the reverse is true.

## 4. NUMERICAL SOLVER: MOLPOP-CEP

In NLTE, the densities of the H₂O states depart from a Boltzmann distribution at the kinetic temperature of the gas. Solving the NLTE problem requires solving the equivalent of Eq. 9 for the densities jointly with Eq. 2 for the radiation field. We have used for that the publicly-available MOLPOP-CEP code (Asensio Ramos & Elitzur 2018), built upon the Coupled Escape Probability (CEP) formalism that tackles exactly the joint population-radiation problems (Elitzur & Asensio Ramos 2006). MOLPOP-CEP outperforms the Λ-iteration technique for short characteristics and parabolic interpolation of the source function, which is a standard NLTE solver for benchmarking. For example, Tables 1-2 and Figures 3-6 in Elitzur & Asensio Ramos (2006) quantify both the computational cost and error incurred by both approaches for a variety of configurations, showing that for a comparable computational cost the error associated with the MOLPOP-CEP calculations is typically orders of magnitude lower. Similar findings have been reported by Yun et al. (2009) in a comparison of the CEP formalism generalized to spherical-shell atmospheres against Monte Carlo simulations.

MOLPOP-CEP requires as input the total density of the target molecule, here o- or p-H₂O. It also requires the list and densities of the colliders for collisions with the target molecule (H₂, H₂O and e⁻ in our application). The model implicitly assumes that the velocities and state-resolved densities of the colliders do not significantly depart from Maxwell-Boltzmann distributions and that the collisional properties for each target-collider pair can be described by a kinetic temperature T. We will relax some of these simplifications in future work to include for example the important effect of non-thermal electrons in the net cooling of the atmosphere (García Muñoz 2023a,b; Gillet et al. 2023). The above properties are specified over a spatial grid of coordinate z that represents the altitude in the atmosphere. As boundary conditions, MOLPOP-CEP requires the external radiation entering the simulation domain. We have considered external irradiation from above and below to mimic the radiation from the star hosting the planet and the radiation upwelling from deeper layers of the atmosphere. We have modified the public version of MOLPOP-CEP by extending the collisional and radiative properties in the model. We have also modified the code to write out the terms that appear in Eqs. 7 and 10, with which we calculate the net cooling rate. The Appendix and Supplementary Information give extended accounts of how MOLPOP-CEP was adapted to our application.

For completeness, we note some potential caveats in the solution of the H₂O NLTE problem with MOLPOP-CEP. Firstly, MOLPOP-CEP assumes a plane-parallel atmosphere. The simplification is likely not critical as we are mainly concerned with a narrow atmospheric layer where the H₂O molecule remains undissociated. Secondly, MOLPOP-CEP describes the radiative transitions of the two H₂O isomers by non-overlapping Gaussian line shapes. The simplification is reasonable at low pressures for which the line wings remain narrow. Thirdly, MOLPOP-CEP cannot accommodate near-resonant channels such as H₂O(010)+O₂(0)↔H₂O(000)+O₂(1) that are important in O₂-rich atmospheres (Yankovsky et al. 2011; Funke et al. 2012; Lang et al. 2020). We assume that such collisional channels are not dominant in the atmospheres that we are exploring, an idea that can always be assessed a posteriori. Lastly, MOLPOP-CEP assumes a static atmosphere. In principle, this simplification may prevent the escape of radiation that occurs when the emitting and absorbing lines become mutually Doppler-shifted in an accelerating gas. We argue in §6.3.7 though that this simplification is unlikely to cause significant errors in the net cooling of the region where the H₂O molecule remains undissociated.
## 5. SOME ESTIMATES

Equation 7 is the basis for our discussion on the net cooling rates. Some insight into it can be gained with a simple two-state molecular model, for which:

$$ 
\begin{aligned}
n_1 + n_2 &= [\text{H}_2\text{O}] \\
\frac{n_2}{n_1} &= \frac{C_{21}}{A_{21}p_{21} + C_{21}g_1} \frac{g_2}{e^{-(E_2-E_1)/kT}} \\
\Gamma^{\text{BB}} &= n_2 A_{21}p_{21}(E_2 - E_1),
\end{aligned}
$$

where we have used detailed balancing to relate the collisional frequencies. Frequent collisions drive the $n_2/n_1$ ratio towards a Boltzmann distribution, with $C_{21}^* \sim A_{21}p_{21}$ setting the condition on the critical density of the background gas and $C_{21} \gg C_{21}^*$ the condition at which LTE is attained. The net radiative bracket $p_{21}$ contains information on two non-local phenomena, namely line opacity and the irradiation of the gas by photons emitted elsewhere. Assuming the first of these dominant, $p_{21}$ is often interpretable as the probability that the line photons may escape from their local environment (Irons 1978), in which case $p_{21} \ll 1$ if the line is optically thick and $\sim 1$ if optically thin. In these limits, it is said that self-absorption is strong and weak, respectively. A small $p_{21}$ helps drive the states towards LTE, possibly ensuring a large $n_2/n_1$ ratio. As a competing effect, a small $p_{21}$ causes the attenuation of the net cooling rate by the effective reduction of the transition probability. These effects combine very differently in purely rotational and in ro-vibrational lines.

Purely rotational transitions between two states with $v=(000)$ are typically characterized by moderate transition probabilities, small energy differences and large densities. For them $p_{21}$ becomes $\ll 1$ in the moderately dense layers of the atmosphere. Ro-vibrational transitions have in comparison large transition probabilities and energy differences, but the density $n_2$ is likely small if the state is populated through collisions and the temperatures remain moderate, i.e. if $(E_2-E_1)/kT \gg 1$. The line opacities associated with ro-vibrational transitions are generally much smaller than those of purely rotational transitions and their net radiative brackets are closer to 1.

### 5.1. Critical densities for rotational and vibrational thermalization

We use the data compiled in Appendix A (the band averages of the transition probabilities and rate coefficients for deexcitation) to estimate some critical densities for the background gas. We do it for an arbitrary $T=400$ K to offer some initial insight. Considering first rotational thermalization within $v=(000)$, with $A_{v'v''} \sim 1.5$ s$^{-1}$, $k_{v'v''}^{e-} \sim 5 \times 10^{-7}$ cm$^3$s$^{-1}$, $k_{v'v''}^{H_2} \sim 2.5 \times 10^{-10}$ cm$^3$s$^{-1}$, $k_{v'v''}^{H_2O} \sim 1.3 \times 10^{-9}$ cm$^3$s$^{-1}$ and an ad-hoc radiative bracket of 1, the critical densities are as low as $6 \times 10^9$ and $1 \times 10^9$ cm$^{-3}$ for neutral atmospheres of H$_2$ and H$_2$O, respectively. For collisions with electrons, the critical density is $3 \times 10^6$ cm$^{-3}$, well below the peak densities predicted in the upper atmospheres of many exoplanets (see Fig. 2 of García Muñoz et al. 2021). In such cases, the collisions with electrons ensure the rotational thermalization of H$_2$O to lower pressures than enabled by the collisions with neutrals alone.

Listed in Table 1 are the critical densities for the vibrational states $v=(010)$, (020), (100) and (001). We have estimated them by means of the generalized condition $\sum_{v''<v'} A_{v'v''} \sim \sum_{v''<v'} C_{v'v''}$. Focusing on $v=(010)$, the critical densities of H$_2$ and H$_2$O are $1.5 \times 10^{13}$ and $4.3 \times 10^{11}$ cm$^{-3}$, respectively. Therefore vibrational thermalization of the bending mode occurs at pressures lower by 1-2 orders of magnitude when the background gas is H$_2$O than when it is H$_2$. The corresponding critical density of electrons is $1.7 \times 10^9$ cm$^{-3}$. This is larger than the electron densities usually predicted for the upper atmospheres of exoplanets. The contribution of the electrons to the vibrational thermalization of H$_2$O in the lower atmosphere may also be minor, as the overall fractional ionization decreases towards the lower atmospheric layers and the thermalization is dominated by the neutrals. In other words, except possibly in extreme ionization conditions, which may not be of interest because the amount of H$_2$O left is small, in general the collisions with electrons will not suffice to drive the H$_2$O molecule to vibrational LTE. The other vibrational states offer similar conclusions. According to Table 1, the critical densities for the stretching modes $v=(100)$ and (001) differ by factors of a few. In reality, they are similar due to rapid interconversion collisions.

The critical densities quoted in Table 1 must be viewed as upper limits. Self-absorption of the photons, especially in H$_2$O-rich atmospheres, will effectively reduce the probability of the upper state in the lines to radiate thereby reducing
the critical densities. Because opacity effects are stronger in purely rotational transitions than in ro-vibrational transitions, the reduction will affect rotational thermalization more severely than vibrational thermalization.

It is worth considering whether the departure of the H₂O molecule from LTE might be detectable with transmission spectroscopy. Focusing again on $v=(010)$, which acts as the main reservoir of vibrationally excited H₂O, the critical densities listed in Table 1 translate into atmospheric pressures of ~1 and 2×10⁻² dyn cm⁻² for the LTE-NLTE transition in H₂-dominated and H₂O-rich atmospheres, respectively. These pressures are well below the pressures of mbars (~10³ dyn cm⁻²) that are often quoted to be probed with low-resolution transmission spectroscopy, which suggests that the detection of H₂O in NLTE with for example JWST is hopeless. The technique of high-resolution spectroscopy (Snellen 2014; López-Morales et al. 2019), which probes the cores of the molecular lines, might be better suited to identify the departure of H₂O from LTE. Dedicated simulations are needed to confirm this possibility.

## 5.2. Opacity in purely rotational and ro-vibrational lines

From Eq. 4, the opacity or optical thickness of a bound-bound radiative transition in a uniform gas is approximately:

$$ \tau_{ij} = \frac{h\lambda_{ij}}{4\pi\Delta\lambda_{ij}}(N_i B_{ij} - N_j B_{ji}), $$

where $\Delta\lambda_{ij}$ is the line width, and $N_i$ and $N_j$ the gas columns along the line of sight (assuming the densities $n_i$ and $n_j$ remain constant). Taking for $\Delta\lambda_{ij}$ the full width at half maximum (fwhm) for Doppler broadening:

$$ \tau_{ij} = \frac{1}{16\pi\sqrt{2\ln 2}} \sqrt{\frac{m}{kT}} A_{ji} \lambda_{ij}^3 (N_i \frac{g_j}{g_i} - N_j), \quad (11) $$

and therefore the line opacity depends to first order on $A_{ji}\lambda_{ij}^3$. The data in Appendix A suggest that for two lines sharing the same lower state, the opacity of the purely rotational line can be orders of magnitude larger than the opacity of the ro-vibrational line. The reverse is expected for the net radiative brackets, which behave approximately as $p_{ji} \sim 1/\tau_{ij}$ for $\tau_{ij} \gg 1$. We have confirmed with the numerical solutions of our NLTE model that indeed the $p_{ji}$ values are typically notably smaller for purely rotational lines than for ro-vibrational ones.

## 5.3. Pressure range to be investigated

Generally, the mass loss of close-in exoplanets is driven by the deposition of stellar extreme-ultraviolet (XUV; wavelengths shortwards of the 912-Å Lyman continuum threshold) photons. For photons of a given wavelength, energy deposition peaks where the optical thickness $\tau_{XUV}=\sigma_{XUV}N_{XUV} \sim 1$, with $\sigma_{XUV}$ being the appropriate cross section and $N_{XUV}$ the overhead column of the XUV-absorbing molecule. This translates into the condition $p_{\tau_{XUV}=1} \sigma_{XUV}/\mu g \sim 1$, which shows that the $\tau_{XUV} \sim 1$ level shifts to high pressures for atmospheres of high molecular weight $\mu$ or under the effect of large gravitational accelerations $g$. For the fiducial values $\sigma_{XUV} \sim 10^{-17}$ cm², $\mu$ in the range 2-18 amu and $g/g_{\oplus}$ in the range 1-5, it is obtained that $p_{\tau_{XUV}=1} = 3.3 \times 10^{-4} - 1.5 \times 10^{-2}$ dyn cm⁻². This defines approximately the region where most of the XUV energy is deposited and where the cooling by H₂O might be more effective at cancelling that effect. Hereafter, we refer to this as the XUV-energy-deposition region.

# 6. SOLUTIONS. ATMOSPHERES THAT ARE NOT EXTERNALLY IRRADIATED

We have run a number of H₂O NLTE simulations of idealized atmospheres to understand the main features of the cooling rate and its dependence on temperature and gas composition. For the simulations described in this section, the atmosphere is not irradiated from the outside and the only radiation that occurs is that arising internally. This condition is relaxed for the simulations in §7. To keep the exercise simple yet informative, we have considered the gas to be isothermal and that its pressure drops exponentially with a scale height $H=kT/\mu g$. The simulation domain (the 'atmosphere') is discretized over a uniform spatial grid of coordinate $z$ such that each scale height is resolved by two bins. The grid spans 16 scale heights from approximately $10^2$ to $10^{-5}$ dyn cm⁻² in 32 bins. We performed a few additional simulations with four bins per scale height. They showed differences in the net cooling rates relative to the standard grid implementation that were typically <1 % and acceptable for our purposes. We have adopted for the gravitational acceleration the terrestrial value $g_{\oplus}=980$ cm s⁻², and explored kinetic temperatures of the background
gas of 200, 800 and 1,500 K. The abundances of H₂, H₂O and e⁻ affect the molecular weight µ, and thus the conversion between z and pressure, as well as the collisional-radiative properties. We have described the gas composition in terms of two parameters (assumed constant throughout the atmosphere), namely: the ratio of densities of the neutrals [H₂O]:[H₂] and the fractional ionization [e⁻]:([H₂]+[H₂O]). We have run separate simulations for o- and p-H₂O using the canonical [o-H₂O]:[p-H₂O]=3:1. The ratio Γº,BB/Γp,BB between their net cooling rates lies typically between ~1 and ~3 due to a combination of abundance and opacity effects. This point is further discussed in §6.3.4. Unless noted otherwise, we focus on the cooling rate Γº,BB due to o-H₂O. We have calculated the cooling rates by means of both Eqs. 7 and 10 using the converged MOLPOP-CEP solutions and confirmed that in the vast majority of simulations they match to within 3-4 significant digits. In the presentation of the cooling rates, we use the units in which they are directly calculated. This choice facilitates the comparison with previous work and is the natural format for implementation in Eq. 1. We note however that in applications concerned with low Mach number flows, as in the study of the terrestrial atmosphere (López-Puertas & Taylor 2001; Feoflov & Kutepov 2012), the cooling rate is often expressed in units of K/day to convey the expected temperature response of the atmosphere under the prescribed forcing.

## 6.1. *Comparison to previous work.*

Neufeld & Melnick (1987) calculated the net cooling rates due to purely rotational transitions within the o-H₂O ground vibrational state for a gas that is not externally irradiated and compared them to similar work by Hollenbach & McKee (1979). We have extracted from our calculations to be discussed in §6.2 those that assume inputs consistent with the conditions in the Neufeld & Melnick (1987) work to prepare their Table 4. Both our calculations and those by Neufeld & Melnick (1987) are presented in Table 2 in the form of cooling efficiencies Γº,BB/[o-H₂O][H₂]. We have adopted the [o-H₂O]:[H₂] ratios indicated there. For each of our entries in the table, the left and right quantities (separated by a comma) are Γº,BB and Γº,BB, respectively.

The cooling efficiencies based on ΓºBB are more relevant for comparison with the Neufeld & Melnick (1987) calculations. Both sets are acceptably consistent although they differ by factors of a few for the higher temperatures and larger η values explored here. As defined in Table 2, η quantifies amongst other things the depth within the gas of the radiating parcel. Also, the level of agreement is sensitive to the adopted [o-H₂O]:[H₂] ratio. It is difficult to identify the reasons for the differences, but they are likely caused by one or more of the following possibilities: 1) Our formulation considers various vibrational states, whereas Neufeld & Melnick (1987) consider only the ground vibrational state; 2) The implemented collisional rate coefficients are different, which should have a stronger impact at low densities when the states do not follow a Boltzmann distribution; 3) Our formulation refers to a gas with an exponential decay in density and zero velocity, whereas Neufeld & Melnick (1987) refer to a uniform-density column with a large velocity gradient; furthermore, it is unclear how accurate is the conversion between the Neufeld & Melnick (1987) formulations for zero and large velocity gradients, which should affect the mapping of the calculations with η.

The comparison between our cooling efficiencies based on Γº,BB and ΓºBB highlights the importance of vibrational excitation and opacity. Indeed, Table 2 shows that all three possibilities Γº,BB>,≈,<ΓºBB occur. The inequality > occurs when the net effect of the vibrationally excited states is cooling the gas. In turn, < occurs when one or more of the ΓBB take negative values, which indicates that some of the ro-vibrational lines absorb radiation rather than emit it. The latter occurs when the gas is heated by radiation emitted elsewhere (typically, but not necessarily, from a deeper layer) yet the gas cannot fully radiate away that energy, which becomes locally deposited. The situation is intimately connected to the non-local nature of the problem, as described further in §6.3.2.

## 6.2. *Calculations*

Figures 2-4 summarize our calculations for atmospheres that are not externally irradiated. Each panel contains multiple curves. Black curves show the net cooling rates Γº,BB. The red symbols indicate ΓºBB with v=(000) (circles), (010) (diamonds), (020) (triangles), (100) (squares), and (001) (stars). Filled and open symbols are used to distinguish rates >0 (cooling) and <0 (heating). For making more explicit the roles of NLTE and opacity, the magenta curves indicate the cooling rates for three idealizations that translate in three variations of Eq. 7. The ‘LTE-thin’ idealization replaces in Eq. 7 the densities obtained from the NLTE solution by the densities in LTE, i.e. $n_j \rightarrow n_j^{LTE}$, and the radiative brackets $p_{ji} \rightarrow 1$; the ‘NLTE-thin’ idealization replaces only $p_{ji} \rightarrow 1$; the ‘LTE-no-thin’ idealization replaces
only $n_j \rightarrow n_j^{\text{LTE}}$.

We begin the discussion with the simulations for a neutral atmosphere with [H₂O]:[H₂]=10⁻⁴ in Fig. 2. In these conditions the H₂O molecule is mostly excited in collisions with H₂ and opacity effects are weaker than for higher H₂O abundances. The calculations for T=200 K are the easiest to understand. Rotational cooling dominates over vibrational cooling and Γ₀,BB ≈ Γ₀,BB₀₀ throughout the atmosphere. For such a low temperature collisions cannot overcome the difference in energy from the ground to the excited vibrational states. The fact that at high pressure Γ₀,BBNLTE-thin ≈ Γ₀,BBLTE-thin ≫ Γ₀,BB reveals the importance of opacity, which attenuates the cooling rate by orders of magnitude with respect to optically thin conditions. For p ≳ 10⁻³ dyn cm⁻² the rotational states within v=(000) are to practical effects thermalized and Γ₀,BB ≈ Γ₀,BBLTE-no-thin. At lower pressures, Γ₀,BB approaches but remains consistently smaller than ΓBBNLTE-thin. The reason is that the net radiative brackets p<sub>ji</sub> do not tend to 1, as might be expected from considerations based on the overhead H₂O column. Rather, p<sub>ji</sub> < 0 for some strong lines, a telltale for upwelling radiation that originates from the deep layers and becomes deposited at low pressures after traversing a sufficiently large gas column. The increased cooling at p ≈ 10² dyn cm⁻² is a general feature of all simulations and results from keeping the bottom boundary of the simulation domain open to outgoing radiation. This feature does not affect the solution at lower pressures and might be avoided by extending the spatial grid to higher pressures than needed.

The calculations for T=800 and 1,500 K reveal additional properties. At high pressures, the cooling is dominated by ro-vibrational lines connecting the ground and first excited vibrational states. The reason is that as the temperature increases the collisional excitation of the vibrational states becomes energetically viable and vibrational cooling efficient. Additionally (see §5.2), the smaller opacities of the ro-vibrational lines facilitate the escape of the emitted photons. In contrast, at low pressures purely rotational cooling is dominant as self-absorption becomes weak for all lines.

Comparison of the three panels for [H₂O]:[H₂]=10⁻⁴ in Fig. 2 shows that the departure of the net cooling rate from the LTE-no-thin idealization occurs at higher pressures for the high-temperature simulations. Equivalently, NLTE affects more the net cooling rate at high temperatures. The trend is a consequence of the increasing contribution of vibrational cooling at high temperatures and the comparative difficulty of the excited vibrational states to thermalize.

The other panels in Figs. 2-3 extend the calculations for neutral atmospheres up to [H₂O]:[H₂]=10², in which limit the atmosphere is virtually pure H₂O. The cooling rates can be for the most part understood with the arguments introduced above, and only a few complementary points are noted here. Firstly, as the H₂O abundance increases the emitted photons become easily absorbed by other H₂O molecules. Increased opacity sets a limit on the net cooling rate, which causes the progressive flattening of the Γ₀,BB curves as [H₂O]:[H₂] increases. Secondly, as the H₂O abundance increases, self-collisions facilitate the excitation of the H₂O states. The latter effect, together with the increased availability of H₂O, dominates over the competing effect of self-absorption and causes that the net cooling rate keeps increasing at even large [H₂O]:[H₂] ratios. Lastly, in the XUV-energy-deposition region, Γ₀,BB attains values ≳ 10⁻⁴ erg cm⁻³s⁻¹ for [H₂O]:[H₂] ≳ 1 and T ≳ 800 K. A significant contribution to it comes from Γ₀,BB₀₁₀.

The extent to which the net cooling rates predicted here can offset the stellar heating that drives the mass loss of close-in exoplanets will depend on the specifics of the planet-star system (e.g. atmospheric composition and stellar spectrum) and must be investigated on a case-by-case basis. Such an investigation is currently under way and will be presented elsewhere. For the time being, we note that XUV heating rates on the order of 10⁻⁶-10⁻⁵ erg cm⁻³s⁻¹ have been reported for hot Jupiters and sub-Neptunes (García Muñoz 2007; García Muñoz et al. 2020), which suggests that H₂O might contribute to the efficient cooling of upper atmospheres provided it survives in sufficient amounts.

The calculations in Fig. 4 extend some of those in Figs. 2-3 by considering a fractional ionization of 10⁻³. In principle, the enhanced collisional excitation caused by the electrons without additional self-absorption should increase the cooling rate. The effect appears however subtle for the conditions explored here. At low temperatures, the purely rotational states remain thermalized up to low pressures without the need for electrons, and vibrational excitation remains energetically unfeasible. At low temperature therefore the electrons modify little the overall population of H₂O states and the cooling. As the temperature increases, rotational cooling remains dominant at low pressures but vibrational cooling takes over at high pressures. The contribution of the electrons to vibrational cooling is modest
because, focusing on $v=(010)$, its excitation rate coefficient in collisions with electrons decreases with temperature whereas the rate coefficients for collisions with H₂ and H₂O either increase with temperature or remain essentially temperature-independent. For example, the ratio of band rate coefficients for (de)excitation of $v=(010)$ in collisions with electrons and H₂ is $k_{v'v''}^{-e}/k_{v'v''}^{H_2} \approx 16,000$ at 200 K but ~540 at 1,600 K (Appendix A). The ratios for collisions with electrons and H₂O are ~380 and ~120. Clearly, a fractional ionization of $10^{-3}$ is insufficient to make a huge impact on the vibrational cooling by H₂O. This finding is likely general in that electrons, under the conditions of planetary atmospheres, contribute to the cooling rate by H₂O but rarely in a dominant way. The idea should be tested with more realistic atmospheres.

The H₂ molecule is an inefficient coolant when it emits at IR wavelengths from its ground electronic state. Using the efficiencies reported by Coppola et al. (2019), we estimate net cooling rates due to H₂ at $p \sim 100$ dyn cm$^{-2}$ of $\sim 10^{-8}$ and $10^{-5}$ erg cm$^{-3}$s$^{-1}$ at temperatures of 200 and 1,500 K, respectively. Because opacity effects remain negligible for the H₂ IR bands, the cooling rate will drop in the same way as pressure for $p<100$ dyn cm$^{-2}$. Compared to our calculations for H₂O in Figs. 2-4, the net cooling rates by H₂ are, even for [H₂O]:[H₂]=10$^{-4}$, lower by at least two orders of magnitude than what our calculations show for H₂O.

## 6.3. Additional insight into the solutions

### 6.3.1. LTE-NLTE transition

It is instructive to quantify the departure of the H₂O states from the Boltzmann distribution that defines LTE. We have produced the ratios $n_i/n_i^{LTE}$ as a function of the energies $E_i$ for three of the simulations discussed above with $T=800$ K. Two simulations are for neutral atmospheres with [H₂O]:[H₂]=10$^{-4}$ and 1, respectively. The third one has a fractional ionization of 10$^{-3}$ and [H₂O]:[H₂]=1. The top three rows of Fig. 5 contain a total of six panels, which correspond to the three simulations at pressures $\approx 10^{-2}$ (left column) and 10$^{-3}$ (right column) dyn cm$^{-2}$.

In all cases the rotational states within $v=(000)$ are the first to thermalize and attain $n_i/n_i^{LTE} \approx 1$, in particular the states with low excitation energy. Rotational thermalization within the other vibrational states is also fast and causes the ratios $n_i/n_i^{LTE}$ to become weakly dependent on $E_i$ yet significantly <1. In other words, the H₂O molecule remains in vibrational NLTE, which translates into the departure of the net cooling rate from the LTE-no-thin idealization in Figs. 2-4. As expected from the combined effects of enhanced collisional frequency and line opacity, a higher pressure or the increased abundance of H₂O or electrons drive the H₂O states towards LTE. The information contained in Fig. 5 complements the discussion on critical densities in section §5.1 and Table 1.

### 6.3.2. Cooling vs. absorption

Some of the calculations in Figs. 2-4 show that the $\Gamma_v^{o, BB}$ become negative (open symbols) at low pressures, and therefore associated with heating, a result that is traced to the net radiative brackets becoming locally negative for some strong lines. The top left panel of Fig. 6 offers further insight into this phenomenon for a neutral atmosphere with [H₂]:[H₂O]=1 and $T=800$ K. Specifically, we have run simulations in which the bottom boundary was progressively shifted to lower pressures. It turns out that $\Gamma_{010}^{o, BB}$ (diamonds in the figure; other $\Gamma_v^{o, BB}$ omitted for clarity) flips from negative to positive when the simulation domain excludes the atmosphere with $p \gtrsim 10^{-3}$ dyn cm$^{-2}$. We argue from this that $\Gamma_v^{o, BB} < 0$ is caused by upwelling radiation that, after traversing a sufficiently large column, becomes deposited in a layer where collisional excitation is inefficient and cannot radiate away the upwelling photons. The panel shows also that although the $\Gamma_v^{o, BB}$ vary between simulations, the net cooling rates seem relatively unaffected by the size of the simulation domain. Additional calculations presented in the top right panel support this interpretation. In their preparation, we kept unchanged the density profiles but reduced the temperature in the lower layers with the goal of diminishing the radiation emitted from them. The results show that when the temperature in the lower layers becomes low enough, $\Gamma_{010}^{o, BB}$ flips from negative into positive.

Related to the above, we have confirmed that the NLTE solutions converge to the NLTE-thin idealization under conditions that minimize the upwelling radiation and its deposition. For a neutral atmosphere with [H₂O]:[H₂]=10$^{-4}$ (and therefore reduced opacity) and $T=800$ K, the bottom panel of Fig. 6 shows that indeed both curves overlap when the geometrical extent of the atmosphere is small and the high-pressure atmospheric layers are omitted (both conditions reducing the upwelling radiation).
### 6.3.3. Strong lines

For reference, we have compiled in Table 3 some information about the strongest lines in a neutral atmosphere at two pressures. The information is organized by (decreasing) line strength or energy exchanged through the line. For $p \approx 0.15$ dyn cm$^{-2}$, most of the lines are cooling and therefore only cooling lines are tabulated. They are dominated by ro-vibrational transitions with upper states in $v=(010)$ or $(020)$. For $p \approx 3 \times 10^{-5}$ dyn cm$^{-2}$, both cooling and heating lines are listed. The cooling lines involve a combination of ro-vibrational transitions with upper states in $v=(010)$ and purely rotational transitions within the ground vibrational state. The heating lines involve ro-vibrational transitions between $v=(000)$ and $(010)$.

### 6.3.4. Isomer-specific cooling

Figure 7 shows the ratio $\Gamma^{o,BB}/\Gamma^{p,BB}$ of the net cooling rates due to o- and p-H₂O for some of the simulations. Typically, the ratio →1 at high pressure and →3 at low pressure, a behavior dictated by line opacity. In the lower layers of the atmosphere, the net radiative bracket behaves overall as the inverse of the density and both contributions largely cancel out in Eq. 7. In such conditions, both isomers are expected to contribute about equally, which is what is found. In contrast, opacity effects are weak at low pressures in the atmosphere and the net cooling rates scale approximately with the relative abundance of the isomer. The transition between the two limits depends on the H₂O abundance and temperature. Higher [H₂O]:[H₂] ratios and lower temperatures enhance the significance of opacity, as discussed above, and this results in the ratio $\Gamma^{o,BB}/\Gamma^{p,BB} \sim 1$ reaching up to lower pressures. It is not evident how to trace the transition between the two limits, which suggests that the calculation of the total cooling rate from H₂O requires separate calculations of the cooling rates from each isomer.

### 6.3.5. Internal energy of excited states

The internal energy of the gas is calculated by summing up the internal energies of all the particles over all the ionization and excitation states within each particle. In NLTE conditions, the H₂O excitation energy contributes by $\rho e_{exc} = \sum_i n_i E_i$ to the internal energy of the gas and to the total energy $\rho E$. We have calculated $\rho e_{exc}$, splitting the summation into purely rotational and ro-vibrational contributions, and compared them to the translational energy $\rho e_{tra} = 3 \sum_i n_i kT/2$. Figure 8 summarizes this comparison for some simulations. It appears that the purely rotational contribution is dominant over the vibrational contribution in the range of temperatures explored, regardless of whether the H₂O molecule is locally in LTE or NLTE. The vibrational energy is by no means negligible though, especially at the higher temperatures and pressures. At $p \gtrsim 0.01-0.1$ dyn cm$^{-2}$ both the rotational and vibrational contributions appear as plateaus characteristic of LTE, which means that the total excitation energy could be estimated from the usual statistical thermodynamics prescriptions (Zel'dovich & Raizer 2002, Chapter III). At lower pressures, the excitation energy should ideally be based on the NLTE solutions. These considerations must be taken into account in the evaluation of the internal energy of the gas when solving Eq. 1.

### 6.3.6. Photodissociation of vibrationally excited states

Photodissociation by FUV (912–2000 Å) photons, with some contribution from XUV photons, determine the destruction of the H₂O molecule in a variety of atmospheres. The existent measurements of H₂O photodissociation cross sections have been obtained at ambient temperature, in which condition most of the H₂O states are in the ground vibrational state. We are unaware of any measurements or theoretical calculations of the photodissociation cross sections at warm temperatures or from the excited vibrational states included in our molecular model. For the atmospheres discussed in §6.3.5, we have calculated the fraction of the H₂O molecules in vibrationally excited states and shown them in Fig. 9. In the NLTE region at low pressures, the population of excited vibrational states is dominated by $v=(010)$, as in the LTE prediction, yet the population is significantly below the LTE prediction. The photodissociation cross section of H₂O in the FUV is on the order of a few times $10^{-18}$ cm² (Hrodmarsson & van Dishoeck 2023). Based on Fig. 9, the photodissociation cross sections of the vibrationally excited states, in particular $v=(010)$, should be about 20-80 times larger than the ground state cross sections to contribute comparably to the H₂O photodissociation at 800 K. The factor drops to 5-14 at 1,500 K. If future calculations or measurements of the photodissociation cross sections from the excited vibrational states prove to be this large, the relevant photodissociation channels should be included in the formulation of the NLTE problem.

### 6.3.7. Static vs. dynamical atmosphere
The acceleration of the gas in the atmosphere Doppler-shifts the emitted photons between their emission and absorption sites, facilitating their escape once the shift becomes comparable to the line width. For a molecular weight $\mu$=18 amu at a fiducial temperature of 400 K, the Doppler fwhm of a H₂O line is ~1 km s⁻¹. Velocities of this order are predicted for the escaping atmospheres of exoplanets, but the key question is whether they are reached in the XUV-energy-deposition region where H₂O contributes the most to the energy balance of the upper atmosphere.

MOLPOP-CEP does not include such a dynamical effect in the treatment of the radiative transfer equation. Consequently, the cooling rates reported here should be viewed as lower limits of what might occur in the atmosphere. Ideally, the impact of this simplification must be assessed on a case-by-case basis. As a first example, we have done such an assessment with the simulations of exoplanet Trappist-1 b presented by García Muñoz (2023b). We have calculated the true overhead column of H₂O at each point of the model spatial grid, $N(H_2O)$. We have also calculated the large-velocity-gradient column $N_{\text{lvg}}(H_2O)$ obtained from equating Eqs. 4 and 5 in Neufeld & Melnick (1987). $N_{\text{lvg}}(H_2O)$ represents the H₂O column in a virtual static atmosphere that would result in the same attenuation of the cooling rate as in the accelerating gas. If $N(H_2O) < N_{\text{lvg}}(H_2O)$ the true H₂O column dictates the line opacity, whereas if $N(H_2O) > N_{\text{lvg}}(H_2O)$ the Doppler shift prevails and the line opacity falls below what is predicted in a static atmosphere. We find that for the profiles of velocity and H₂O density predicted in the Trappist-1 b atmosphere, $N(H_2O) \lesssim N_{\text{lvg}}(H_2O)$ throughout the region where H₂O remains undissociated. In this specific example the Doppler-shifting of the lines should not significantly affect the escape of the photons and the net cooling rate with respect to the calculations for a static atmosphere.

## 7. SOLUTIONS. ATMOSPHERES THAT ARE EXTERNALLY IRRADIATED

A real atmosphere is irradiated from below by photons emerging from deeper layers and from above by stellar photons. For bound-bound radiative transitions, the molecule is initially photoexcited to an upper state. Table 4 (Appendix A; last two columns) lists the rates at which the energy is pumped from $v$=(000) under optically thin conditions. They depend on the convolved effect of the H₂O radiative properties and the stellar spectrum, and we show them for two stellar types. In both cases, most of the stellar energy that is pumped from $v$=(000) ends up into $v$=(001). Some of the excited molecules will spontaneously radiate away their energy, not heating the local gas. Others will transfer their energy to the background gas through collisions, heating it. Photoexcitation sets a cascade (Fig. 1) in which the initial photon energy is communicated to the various H₂O states through collisions and radiation. In this section, we consider how external irradiation affects the net cooling rates of some of the atmospheres described above.

MOLPOP-CEP can optionally consider at each boundary of the simulation domain isotropic radiation of radiance $WB_{\lambda}(T_{bb})$ entering the domain. $W$ is the so-called dilution factor and $B_{\lambda}(T_{bb})$ the blackbody radiance at temperature $T_{bb}$. The flux of energy through the boundary into the simulation domain is $\pi WB_{\lambda}(T_{bb})$ (units of erg cm⁻²s⁻¹). At the bottom boundary, we have prescribed $W$=1 and a blackbody temperature that coincides with that of the gas, i.e. $T_{bb}=T$. The treatment is a valid representation of an atmosphere with equilibrium temperature $T_{eq} \approx T$. For completeness, we note that MOLPOP-CEP includes the option of prescribing the external irradiation through a user-prepared table, which might be important at short wavelengths for which the stellar spectrum notably departs from that of a black body. At the top boundary, we have prescribed $W=(R_{\star}/a)^2$, where $R_{\star}$ and $a$ are the stellar radius and planet-star orbital distance, respectively, and for $T_{bb}$ a value of 2,500 K. These choices ensure that the flux of energy through the top boundary is a fair representation of what the real atmosphere receives at IR wavelengths from a TRAPPIST-1-like star that subtends the solid angle $\pi(R_{\star}/a)^2$. To specify the orbital distance, we have resorted to the definition of $T_{eq}$ for a zero-albedo atmosphere without energy redistribution, for which $(R_{\star}/a)^2=2(T_{eq}/T_{\star})^4$ (Traub & Oppenheimer 2010). Therefore $W$ depends only on the ratio between the planet equilibrium temperature and the stellar effective temperature, for which we have adopted $T_{eq}=T$ and $T_{\star}=2,500$ K. For reference, at the top boundary $W=8.2 \times 10^{-5}$, 0.021 and 0.26 for $T_{eq}=200, 800$ and $1,500$ K. Figure 10 sketches some of these ideas.

Figure 11 summarizes the net cooling rates for neutral atmospheres irradiated from above only or simultaneously from above and below. Both sets of curves are represented with the same combination of lines (solid for net cooling; dotted for net heating) and symbols. They are easy to tell apart because irradiation from below reduces the net cooling at $p \gtrsim 1-10$ dyn cm⁻² but leaves the lower pressures largely unaffected. The figure shows for comparison the net cooling
rates for atmospheres that are not externally irradiated (solid line, no symbols). For the conditions explored here, irradiation from above seems to have a larger effect on the net cooling rate than irradiation from below. The region in the atmosphere affected by stellar IR irradiation depends on the H₂O abundance, but it appears that the condition $p \gtrsim 0.1$ dyn cm⁻² serves as a useful reference. This pressure level is largely dictated by the condition that the excited vibrational states transfer their energy to the background gas through collisions, and it is therefore connected to the critical densities discussed in section 5.1 and summarized in Table 1. The higher the ratio [H₂O]:[H₂], the lower the pressure at which the energy is deposited, a consequence of higher opacities and deexcitation collisional frequencies. In general, stellar IR irradiation may reduce the net cooling rate by a factor of ~2 at $p \gtrsim 1$ dyn cm⁻² with respect to the solution for atmospheres that are not externally irradiated, or even cause the net heating of the atmosphere in those atmospheric layers.

The last row of Fig. 5 shows the ratios $n_i/n^{LTE_i}$ for a neutral atmosphere with [H₂O]:[H₂]=1 that is externally irradiated. The main difference with the calculations presented in the same figure for atmospheres that are not irradiated is that the states $v=(100)$ and (001) are overpopulated with respect to LTE as a result of stellar-driven photoexcitation. As expected, increasing pressure drives the population towards LTE.

## 8. SUMMARY AND OUTLOOK

We have built a model to investigate the effects of NLTE in the H₂O-rich atmospheres of exoplanets. The model combines up-to-date information on the collisional and radiative properties of H₂O with a state-of-the-art solver for the NLTE problem. Using idealized atmospheres of constant temperature and pressure dropping exponentially, we have explored how the temperature and relative abundances of H₂, H₂O and e⁻ affect the cooling through bound-bound radiative transitions of H₂O and the LTE-NLTE transition in the atmosphere. The net cooling rate typically increases with increasing H₂O abundance and, more weakly, with increasing fractional ionization. These factors also facilitate the thermalization of the H₂O states with respect to what occurs in a neutral H₂-dominated atmosphere. Water remains in rotational LTE up to the pressures where the stellar XUV energy is deposited and the atmospheric gas is accelerated to space. Attaining vibrational LTE requires higher densities than rotational LTE, on the order of ~10¹³ cm⁻³ for $v=(010)$ in a neutral H₂-dominated atmosphere but a few times ~10¹¹ cm⁻³ in a neutral H₂O-rich atmosphere. We have estimated the critical densities for the LTE-NLTE transition of other vibrational states, as the information might be useful for future searches of LTE departures with transmission spectroscopy. High temperatures facilitate the collisional excitation of the H₂O states and therefore enhance the cooling rate. As the temperature increases, the vibrational states, typically in NLTE, contribute more significantly to the cooling. We emphasize the importance of H₂O-H₂O self-collisions and line opacity in the NLTE treatment. Our simulations that include stellar IR irradiation of the atmosphere show that the H₂O populations are notably different with respect to the case of no external irradiation, especially for $v=(100)$ and (001). Most of the heating caused by stellar IR photons occurs at pressures $\gtrsim 0.1$ dyn cm⁻². This is deeper in the atmosphere than where the stellar XUV photons are deposited, but the exact region may depend on other conditions that need to be explored on a case-specific basis including the stellar type and atmospheric properties. Looking ahead, the net cooling rates reported in this work prompt detailed investigations of the mass loss of H₂O-rich atmospheres including H₂O NLTE. The model presented here, together with the corresponding input and output files, are available from the corresponding author upon reasonable request. The collisional and radiative data will also be made available through the EMAA database https://dx.doi.org/10.17178/EMAA.

## ACKNOWLEDGEMENTS

AAR acknowledges support from the Agencia Estatal de Investigación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) and the European Regional Development Fund (ERDF) through project PID2022-136563NB-I0.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Sketch of our molecular model. The quoted energies correspond to the lowest rotational state (including both isomers) within the vibrational state. For each upper vibrational state, the dotted-line arrows indicate the band with the largest transition probability (big dots) and with the second largest transition probability (small dots) (see $A_{v'v''}$ for $v' \neq v''$, Appendix A). For the lower vibrational state $v=(000)$, the squared-line arrows indicate the photoexcitation bands that carry the most of stellar energy (large squares) and the second most (small squares) (see $g_{v''v}^{\odot}hc/\lambda_{v'v''}$ or $g_{v''v}^{\text{tr1}}hc/\lambda_{v'v''}$ in Table 4). The exchange of energy through collisions (not marked in the sketch) occurs mainly from $v=(000)$ into $(010)$ with some contribution from $v=(000)$ into $(020)$. Our treatment assumes that the total density of H$_2$O is known at each location in the spatial grid $z$ and time $t$. In practice, [H$_2$O] would be determined from solving the continuity equation of the ensemble of H$_2$O states without considering a separate continuity equation for each of the individual states. The outputs of the NLTE model include the net cooling rate and the contribution of H$_2$O to the internal excitation energy.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Cooling rates for atmospheres without external irradiation. Conditions in panel headers. Filled and open symbols refer to rates >0 (cooling) and <0 (heating), respectively. See text for explanation of each curve.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Cooling rates for atmospheres without external irradiation. Conditions in panel headers. Filled and open symbols refer to rates >0 (cooling) and <0 (heating), respectively. See text for explanation of each curve.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Cooling rates for atmospheres without external irradiation. Conditions in panel headers. Filled and open symbols refer to rates >0 (cooling) and <0 (heating), respectively. See text for explanation of each curve.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Departures from LTE for the H₂O states at T=800 K. Left and right columns are for *p*≈10⁻² and 10⁻³ dyn cm⁻², respectively. Top three rows are for atmospheres that are not externally irradiated; bottom row is for atmospheres that are externally irradiated. Conditions regarding fractional ionization and water abundance are specified in each panel.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Top left. Net cooling rate $\Gamma^{o,BB}$ (solid line) and partial component $\Gamma^{o,BB}_{010}$ (diamonds; filled and open symbols for local cooling and heating, respectively) for a neutral atmosphere with [H₂O]:[H₂]=1 and T=800 K without external irradiation. Colors are used to differentiate some of the calculations performed over small simulation domains. Top right. Partial component of the cooling rate $\Gamma^{o,BB}_{010}$ (filled and open symbols for local cooling and heating, respectively) for a neutral atmosphere with [H₂O]:[H₂]=1 without external irradiation. The NLTE calculations assume the same profiles of H₂O and H₂ vs. z as for the calculations of the top-left panel. For the temperature, however, we assumed a discontinuous two-step profile with temperature values as indicated in the figure. Colors are used to differentiate the calculations for each temperature profile. The calculations were done with the molecular model truncated to 200 states. Bottom. Net cooling rate (solid black curve, thin line) for a neutral atmosphere with [H₂O]:[H₂]=10⁻⁴ and T=800 K without external irradiation, as described in section 6.2, and corresponding cooling rate for NLTE-thin idealization (dashed magenta curve, thin line). Each of the black/magenta sets represented by thick lines refer to atmospheres with the same characteristics but geometrically truncated.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Ratio of net cooling rates for the two isomer forms of H₂O. Calculations based on the simulations of §6 for atmospheres that are not externally irradiated.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Rotational, meaning within $v=(000)$, and ro-vibrational components of the excitation energy relative to the translational energy. Calculations based on neutral atmospheres with [H₂O]:[H₂]=1 that are not externally irradiated. Exceptionally, the small symbols are for an atmosphere of T=800 K that is externally irradiated as described in section §7.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Fraction of excited vibrational states relative to the H₂O total density. Same atmospheric conditions as for Fig. 8.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Simulation with MOLPOP-CEP of externally irradiated atmospheres. The radiation from below is implemented as isotropic blackbody radiation upwelling from all directions (solid angle $2\pi$) at the planet equilibrium temperature $T_{\text{eq}}$. The radiation from above is implemented as isotropic blackbody radiation entering from all directions (solid angle $2\pi$). The intensity of this radiation is scaled by a dilution $W$ to ensure that the energy entering the atmosphere from above matches the energy that is received by a planet of equilibrium temperature $T_{\text{eq}}$ orbiting a star of effective temperature $T_{\star}$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Cooling rates for atmospheres with external irradiation (as irradiated by a TRAPPIST 1-like star from above and a black body at the stated temperatures from below; see text). Conditions in panel headers. Solid and dotted lines refer to rates >0 (cooling) and <0 (heating), respectively.</figcaption></figure>
<table><caption>**Table 1.** Critical densities [cm⁻³] at T=400 K for the vibrational thermalization of the quoted states in collisions with one of the H₂, H₂O or electron particles. Here and throughout the manuscript X(Y)=X×10⁷.</caption><thead><tr><th rowspan="2">Collider</th><th colspan="4">Upper vibrational state v'</th></tr><tr><th>(010)</th><th>(020)</th><th>(100)</th><th>(001)</th></tr></thead><tbody><tr><td>H₂</td><td>1.5(+13)</td><td>1.4(+13)</td><td>9.5(+12)</td><td>2.3(+12)</td></tr><tr><td>H₂O</td><td>4.3(+11)</td><td>4.3(+11)</td><td>3.7(+11)</td><td>7.2(+10)</td></tr><tr><td>e⁻</td><td>1.7(+09)</td><td>1.7(+09)</td><td>3.1(+08)</td><td>2.0(+09)</td></tr></tbody></table>

<table><caption>**Table 2.** Cooling efficiency Γ°,BB/[o-H₂O][H₂] [erg cm³s⁻¹] as a function of η (≡N(o-H₂O)[H₂](T/1000)⁻⁰.⁵×1.35×10⁻²⁰, cgs units; N is the column to the top of the atmosphere) and temperature. For each η, the top entry (in bold) refers to Table 4 of Neufeld & Melnick (1987). The pairs of entries immediately below are our calculations, arranged by the adopted [o-H₂O]:[H₂] ratio. Left (regular font) and right (bold) entries are based on the net cooling rate Γ°,BB and the partial component Γ°,BB°°°°, respectively. Limiting η≤10² helps approximately realize the condition established by Neufeld & Melnick (1987) of [H₂]<10¹⁰ cm⁻³ in their calculations.</caption><thead><tr><th rowspan="3">η</th><th rowspan="3">[o-H<sub>2</sub>O]:[H<sub>2</sub>]</th><th colspan="6">Temperature [K]</th></tr><tr><th colspan="6">200</th></tr><tr><th>400</th><th>600</th><th>800</th><th>1000</th><th>1500</th></tr></thead><tbody><tr><td rowspan="6">10<sup>0</sup></td><td></td><td><b>1.3</b>(-24)</td><td><b>4.1</b>(-24)</td><td><b>8.9</b>(-24)</td><td><b>1.6</b>(-23)</td><td><b>2.6</b>(-23)</td><td><b>5.8</b>(-23)</td></tr><tr><td>0.75(-5)</td><td>8.4(-25),<b>8.4</b>(-25)</td><td>1.9(-24),<b>2.1</b>(-24)</td><td>2.7(-24),<b>4.2</b>(-24)</td><td>3.6(-24),<b>7.6</b>(-24)</td><td>4.4(-24),<b>1.1</b>(-23)</td><td>6.9(-24),<b>1.4</b>(-23)</td></tr><tr><td>0.75(-4)</td><td>1.0(-24),<b>1.0</b>(-24)</td><td>2.2(-24),<b>2.7</b>(-24)</td><td>3.1(-24),<b>7.7</b>(-24)</td><td>3.9(-24),<b>2.0</b>(-23)</td><td>4.7(-24),<b>3.9</b>(-23)</td><td>6.7(-24),<b>8.5</b>(-23)</td></tr><tr><td>0.75(-3)</td><td>1.2(-24),<b>1.2</b>(-24)</td><td>2.6(-24),<b>3.2</b>(-24)</td><td>3.7(-24),<b>1.1</b>(-23)</td><td>4.7(-24),<b>3.6</b>(-23)</td><td>5.6(-24),<b>8.6</b>(-23)</td><td>7.7(-24),<b>3.0</b>(-22)</td></tr><tr><td>0.75(-2)</td><td>1.3(-24),<b>1.3</b>(-24)</td><td>3.2(-24),<b>3.7</b>(-24)</td><td>4.7(-24),<b>1.1</b>(-23)</td><td>6.1(-24),<b>3.6</b>(-23)</td><td>7.4(-24),<b>9.2</b>(-23)</td><td>1.1(-23),<b>3.9</b>(-22)</td></tr><tr><td></td><td><b>1.1</b>(-24)</td><td><b>3.8</b>(-24)</td><td><b>8.4</b>(-24)</td><td><b>1.6</b>(-23)</td><td><b>2.5</b>(-23)</td><td><b>5.7</b>(-23)</td></tr><tr><td rowspan="6">10<sup>1</sup></td><td></td><td><b>1.1</b>(-24)</td><td><b>3.8</b>(-24)</td><td><b>8.4</b>(-24)</td><td><b>1.6</b>(-23)</td><td><b>2.5</b>(-23)</td><td><b>5.7</b>(-23)</td></tr><tr><td>0.75(-5)</td><td>6.3(-25),<b>6.3</b>(-25)</td><td>1.6(-24),<b>1.6</b>(-24)</td><td>2.4(-24),<b>2.9</b>(-24)</td><td>3.3(-24),<b>4.5</b>(-24)</td><td>4.1(-24),<b>5.9</b>(-24)</td><td>6.6(-24),<b>8.0</b>(-24)</td></tr><tr><td>0.75(-4)</td><td>8.8(-25),<b>8.8</b>(-25)</td><td>2.0(-24),<b>2.2</b>(-24)</td><td>2.9(-24),<b>4.6</b>(-24)</td><td>3.8(-24),<b>9.5</b>(-24)</td><td>4.5(-24),<b>1.6</b>(-23)</td><td>6.6(-24),<b>3.3</b>(-23)</td></tr><tr><td>0.75(-3)</td><td>1.1(-24),<b>1.1</b>(-24)</td><td>2.5(-24),<b>2.9</b>(-24)</td><td>3.6(-24),<b>8.0</b>(-24)</td><td>4.6(-24),<b>2.2</b>(-23)</td><td>5.5(-24),<b>4.5</b>(-23)</td><td>7.6(-24),<b>1.2</b>(-22)</td></tr><tr><td>0.75(-2)</td><td>1.3(-24),<b>1.3</b>(-24)</td><td>3.2(-24),<b>3.6</b>(-24)</td><td>4.7(-24),<b>1.1</b>(-23)</td><td>6.1(-24),<b>3.4</b>(-23)</td><td>7.3(-24),<b>8.3</b>(-23)</td><td>1.1(-23),<b>3.2</b>(-22)</td></tr><tr><td></td><td><b>8.3</b>(-25)</td><td><b>3.2</b>(-24)</td><td><b>7.7</b>(-24)</td><td><b>1.5</b>(-23)</td><td><b>2.4</b>(-23)</td><td><b>5.5</b>(-23)</td></tr><tr><td rowspan="6">10<sup>2</sup></td><td></td><td><b>8.3</b>(-25)</td><td><b>3.2</b>(-24)</td><td><b>7.7</b>(-24)</td><td><b>1.5</b>(-23)</td><td><b>2.4</b>(-23)</td><td><b>5.5</b>(-23)</td></tr><tr><td>0.75(-5)</td><td>3.9(-25),<b>3.9</b>(-25)</td><td>1.2(-24),<b>1.2</b>(-24)</td><td>2.0(-24),<b>2.1</b>(-24)</td><td>2.8(-24),<b>3.0</b>(-24)</td><td>3.6(-24),<b>3.9</b>(-24)</td><td>6.0(-24),<b>5.6</b>(-24)</td></tr><tr><td>0.75(-4)</td><td>6.4(-25),<b>6.4</b>(-25)</td><td>1.6(-24),<b>1.7</b>(-24)</td><td>2.6(-24),<b>3.1</b>(-24)</td><td>3.4(-24),<b>5.2</b>(-24)</td><td>4.2(-24),<b>7.6</b>(-24)</td><td>6.3(-24),<b>1.4</b>(-23)</td></tr><tr><td>0.75(-3)</td><td>9.2(-25),<b>9.3</b>(-25)</td><td>2.2(-24),<b>2.4</b>(-24)</td><td>3.3(-24),<b>4.9</b>(-24)</td><td>4.3(-24),<b>1.0</b>(-23)</td><td>5.2(-24),<b>1.9</b>(-23)</td><td>7.5(-24),<b>4.7</b>(-23)</td></tr><tr><td>0.75(-2)</td><td>1.2(-24),<b>1.2</b>(-24)</td><td>2.9(-24),<b>3.3</b>(-24)</td><td>4.3(-24),<b>7.9</b>(-24)</td><td>5.7(-24),<b>2.0</b>(-23)</td><td>7.0(-24),<b>4.3</b>(-23)</td><td>1.0(-23),<b>1.3</b>(-22)</td></tr></tbody></table>

<table><caption>**Table 3.** Strongest lines in a neutral atmosphere with [H₂O]:[H₂]=1 and T=800 K. I: Wavelength [µm]; II: Energy rate [erg cm⁻³s⁻¹]; III: Transition: Upper-lower states, in the form (v₁v₂v₃;Jτ). For heating, the energy rate of the line is negative.</caption><thead><tr><th rowspan="3">#</th><th colspan="3">p≈0.15 dyn cm<sup>-2</sup></th><th colspan="3">p≈5×10<sup>-5</sup> dyn cm<sup>-2</sup></th></tr><tr><th colspan="3">Cooling lines</th><th colspan="3">Cooling lines</th></tr><tr><th>I</th><th>II</th><th>III</th><th>I</th><th>II</th><th>III</th></tr></thead><tbody><tr><td>1</td><td>5.24</td><td>7.87(-6)</td><td>(010;18-17)-(000;17-17)</td><td>7.04</td><td>5.04(-8)</td><td>(010;33)-(000;43)</td></tr><tr><td>2</td><td>5.44</td><td>7.81(-6)</td><td>(020;43)-(010;33)</td><td>7.17</td><td>4.71(-8)</td><td>(010;41)-(000;53)</td></tr><tr><td>3</td><td>4.73</td><td>7.81(-6)</td><td>(010;95)-(000;83)</td><td>7.21</td><td>4.52(-8)</td><td>(010;43)-(000;55)</td></tr><tr><td>4</td><td>4.79</td><td>7.81(-6)</td><td>(010;101)-(000;91)</td><td>6.86</td><td>4.33(-8)</td><td>(010;21)-(000;33)</td></tr><tr><td>5</td><td>5.11</td><td>7.79(-6)</td><td>(010;15-11)-(000;14-9)</td><td>6.43</td><td>4.24(-8)</td><td>(010;5-3)-(000;5-1)</td></tr><tr><td>6</td><td>4.68</td><td>7.78(-6)</td><td>(010;87)-(000;77)</td><td>7.33</td><td>4.04(-8)</td><td>(010;5-3)-(000;6-3)</td></tr><tr><td>7</td><td>5.17</td><td>7.78(-6)</td><td>(010;16-13)-(000;15-13)</td><td>35.43</td><td>3.79(-8)</td><td>(000;71)-(000;6-1)</td></tr><tr><td>8</td><td>5.57</td><td>7.78(-6)</td><td>(020;13-13)-(010;12-11)</td><td>32.99</td><td>3.68(-8)</td><td>(000;7-3)-(000;6-1)</td></tr><tr><td>9</td><td>5.36</td><td>7.76(-6)</td><td>(020;5-3)-(010;41)</td><td>30.87</td><td>3.67(-8)</td><td>(000;8-1)-(000;7-1)</td></tr><tr><td>10</td><td>5.62</td><td>7.75(-6)</td><td>(020;12-11)-(010;11-11)</td><td>35.94</td><td>3.62(-8)</td><td>(000;6-3)-(000;5-3)</td></tr></tbody></table>
# APPENDIX

## A. MOLECULAR MODEL

### A.1. *Energies, statistical weights and transition probabilities*

We adopted the energies, statistical weights and transition probabilities for H₂O from the BT2 linelist (Barber et al. 2006), which refers to ~28,500 ro-vibrational states with identified quantum numbers. Because the cost of a NLTE calculation rises steeply with the number of states and convergence difficulties arise when negligibly populated states are present, we made an effort to simplify the molecular model without compromising the accuracy of the NLTE calculations, as described below.

Far ultraviolet (FUV) photons excite the H₂O molecule into various electronic states that tend to predissociate (Smith et al. 1981). These photodissociation channels are generally included in the hydrodynamical models that govern the H₂O density in planetary atmospheres. To avoid double-counting, we omitted from our molecular model the explicit treatment of excited electronic states.

Focusing on the ground electronic state, after some experimentation we settled on a molecular model that includes every ro-vibrational state with quantum numbers $v=(000), (010), (020), (100), (001), J \le 20$ and energies $E_i/k \le 7,200$ K ($k$ is Boltzmann's constant). It is nearly identical to that used by Faure & Josselin (2008) to describe the collisions of H₂O with $e^-$ and H₂, and includes 408 and 410 states of o- and p-H₂O respectively. The molecular model accommodates the majority of the radiative transitions and collisional channels that are significant at the temperatures and irradiation conditions of interest here. For reference, it includes ~160 states of each isomer with up to $J=20$ in $v=(000)$ and ~40 states with up to $J=10$ in $v=(001)$. A more comprehensive molecular model would include about the same number of ro-vibrational states in $v=(001)$ as in (000) to ensure their complete connection by stellar IR radiation. The information in Table 4 confirms that the vibrational states $v=(010), (020), (100)$ and (001) encompass $\ge 90\%$ of the energy pumped from $v=(000)$ into all the excited vibrational states for stellar radiation of effective temperatures 2,500 and 5,780 K. In particular, the energy rates associated with photoexcitation in the optically thin limit are shown as $g_{v''v'}hc/\lambda_{v'v''}$, where $g_{v''v'}$ is the photoexcitation frequency.

The molecular model includes no isotopologues other than H₂¹⁶O. The next most abundant isotopologue in the terrestrial atmosphere is H₂¹⁸O and comprises 0.2% of the total water (Galewsky et al. 2016). Our calculations based on H₂¹⁶O (with the proper scaling of abundances) suggest that for the H₂O columns that we are exploring, the cooling due to H₂¹⁸O may be significant at high pressure, where the overhead column of H₂O is large, but subdominant, especially at high temperatures when vibrational cooling becomes effective.

### A.2. *Collisional (de)excitation.*

We are interested in H₂-H₂O atmospheres that, as a result of dissociation and ionization, may also contain electrons and other heavy particles. These particles undergo inelastic collisions with the target o- and p-H₂O, exchanging energy with them. We compiled up-to-date sets of state-to-state rate coefficients for collisions of H₂O with electrons, H₂, H₂O, H, O, O₂ and He. The NLTE calculations presented here focus on atmospheres composed of H₂, H₂O and $e^-$ only. Compiling the rate coefficients for collisions with the other heavy particles helps anticipate their significance in more complex models. What is meant here by electrons is the population of thermal electrons, assumed to have the same kinetic temperature as the heavy particles. We defer to future work the implementation in the NLTE calculations of non-thermal electrons. We discuss below the rate coefficients for deexcitation. MOLPOP-CEP determines those for excitation internally from application of detailed balancing.

Setting up the NLTE model requires prescribing state-to-state rate coefficients $k_{ji}^Q$ for each deexcitation channel $j \to i$ between individual ro-vibrational states of the o- and p-H₂O target and each projectile or collider Q. This information is not always known, in which cases it must be estimated. We distinguish:

1. When no quantitative information on the rate coefficients exists. In this case, we obtained crude estimates by scaling the rate coefficients for collisions with H₂ as in $k_{ji}^Q=k_{ji}^{H_2}(\mu_{H_2O-H_2}/\mu_{H_2O-Q})^{1/2}$, where $\mu$ (with a
subscript, not to mistake for the molecular weight of the atmosphere) stands for the reduced mass of the collision pair (Walker et al. 2014).

2. When the rate coefficients for rotational deexcitation are known for some states within the ground vibrational state but not for others of typically higher energy. Neufeld & Melnick (1987) devised a method (hereafter, the NM87 method) to estimate the unknown rate coefficients from the known ones. It relies on the partly testable idea that the total relaxation rate coefficient is often weakly dependent on the upper state *j* and may be estimated from the available rate coefficients, i.e. that $k_j^Q = \sum_{i<j} k_{ji}^Q \approx k_u^Q(\downarrow)$ for states *j* above a threshold. We adopted the NM87 method, and for the state-to-state partitioning implemented their $\gamma=0$-scheme (their Eq. 6) so that $k_{ji}^Q/k_u^Q(\downarrow) = g_i / \sum_{m<j} g_m$. We estimated $k_u^Q(\downarrow)$ by averaging the total relaxation rate coefficients $k_j^Q$ over the ~5–10 higher-energy states *j* for which the original data report the rate coefficients. We mostly used the NM87 method to extrapolate, i.e. to estimate the rate coefficients for states *j* above the threshold used in the determination of $k_u^Q(\downarrow)$ but occasionally, in particular for H₂O self-collisions, we used the method also for interpolation when the relevant rate coefficients were missing. In our experience, providing relaxation channels for all the states in the molecular model improves the numerical robustness of MOLPOP-CEP.

3. When the band rate coefficient for deexcitation between vibrational states $k_{v'v''}^Q$ is known, but the state-to-state rate coefficients $k_{ji}^Q$ are unknown. In this case, it is practical to assume that the state-to-state partitioning occurs following the propensity rules of collisions with another heavy particle. As a rule, we interpreted the band rate coefficients from experiments or from calculations that omit rotational details to represent thermal averages of the state-to-state rate coefficients. Formally, the thermal averages are defined as $A_{v'v''}$ in Table 4, replacing $A_{ji}$ by $k_{ji}^Q$, and extending the summations over the relevant states of the molecular model. We further assumed that the branching ratios $t_{ji}^Q = k_{ji}^Q/k_{v'v''}^Q$ are identical to the $t_{ji}^{H_2}$ for collisions with H₂, which we determined from Daniel et al. (2011) and Faure & Josselin (2008) (see below). Given $k_{v'v''}^Q$ and $t_{ji}^Q$, the state-to-state rate coefficients $k_{ji}^Q$ are fully determined. By construction, the thermal averages of $k_{ji}^Q$ are equal to the prescribed $k_{v'v''}^Q$.

Tables 5–7 summarize the band rate coefficients for collisions with electrons, H₂ and H₂O. Some points worth noting from these tables and Table 4 are: 1) H₂O radiates faster when the transition involves a change in vibrational state than when the transition occurs within the same vibrational state; 2) Rotational relaxation in collisions occurs notably faster than vibrational relaxation. For example, in collisions with electrons, relaxation within the ground vibrational state proceeds with a rate coefficient ~5×10⁻⁷ cm³s⁻<fnref n="1" />, whereas the rate coefficient for vibrational relaxation (010)→(000) is ~10⁻⁸ cm³s⁻<fnref n="1" />. Based on this, the first excited vibrational state attains LTE at pressures 2-3 orders of magnitude higher than required for thermalization of the rotational population; 3) Both rotational and vibrational relaxation of H₂O occurs 2-4 orders of magnitude faster in collisions with electrons than with heavy particles. As a consequence, H₂O relaxation is affected by electrons once the gas fractional ionization reaches ~10⁻⁴-10⁻², depending on temperature; 4) Both rotational and vibrational relaxation of H₂O occur more than an order of magnitude faster through self-collisions with other H₂O molecules than through collisions with other heavy particles. Foreseeably, H₂O remains in LTE up to lower pressures in a H₂O-rich atmosphere than in a H₂-dominated atmosphere; 5) The rate coefficients for collisional relaxation depend moderately on temperature. Typically, the relaxation proceeds faster at lower temperatures for collisions with thermal electrons, and at higher temperatures for collisions with H₂ or H₂O.

### A.2.1. Collisions with electrons.

Faure & Josselin (2008) assembled a set of state-to-state rate coefficients for H₂O states in collisions with thermal electrons.<fnref n="1" /> We generally adopted their rate coefficients, as downloaded from the CDS web. For channels (010)→(000), (020)→(000), (020)→(010), (100)→(000), and (001)→(000), we rescaled the data by temperature-dependent factors to ensure that the rescaled thermally-averaged rate coefficients match those reported by Ayouz et al. (2021). The Faure & Josselin (2008) and Ayouz et al. (2021) rate coefficients are generally consistent to within a factor of 3 except for (020)→(000), in which case they differ by more than a factor of 10. The band rate coefficients, calculated as thermal averages of the state-to-state rate coefficients, are shown in Table 5.

<footnote n="1">¹ The data set is currently being upgraded with new rate coefficients for collisions with electrons and H₂ (A. Faure, private communication).</footnote>
### A.2.2. Collisions with H₂.

Phillips et al. (1996), Faure et al. (2006, 2007), Dubernet et al. (2006, 2009), Daniel et al. (2011) and Żóltowski et al. (2021) calculated the rate coefficients for rotational deexcitation within the ground vibrational state of H₂O in collisions with H₂ molecules. We adopted the Daniel et al. (2011) data for collisions with o- and p-H₂ as downloaded from LAMDA (Schöier et al. 2005) and produced rate coefficients relevant to 3:1-ortho:para Boltzmann distributions of H₂. The Daniel et al. (2011) calculations apply to the 45 lower-energy states of o- and p-H₂O. For channels involving higher-energy states, we adopted the rate coefficients reported by Faure & Josselin (2008).

Much less is known about the rate coefficients for channels involving excited vibrational states. The only relevant measurements that exist for collisions of H₂O with H₂ seem to be those by Zittel & Masturzo (1991). Faure & Josselin (2008) summarized the experimental and theoretical information available at the time and produced a set of state-to-state rate coefficients for collisions involving various vibrational states. We adopted their deexcitation rate coefficients except for channels within the ground vibrational state, as described above. Recent calculations for (010)→(000) include those by Stoecklin et al. (2019), Wiesenfeld (2021, 2022) and García-Vázquez et al. (2023), but they do not cover the temperature range investigated here. In particular, the García-Vázquez et al. (2023) work (their Fig. 9) shows that (010) relaxation occurs with a rate coefficient ~1.3×10⁻¹² cm³ s⁻¹ at 300 K that is in agreement with the experiments. Their calculations also support the usual assumption that rotational relaxation occurs rapidly within each vibrational state.

Of importance for describing the H₂O collisions with other heavy particles, we used the combined Daniel et al. (2011) and Faure & Josselin (2008) rate coefficients to define branching ratios $t_{ji}^{H_2} \equiv k_{ji}^{H_2} / k_{v'v''}^{H_2}$, where the numerator refers to the state-to-state rate coefficient between states j and i (in vibrational states v' and v'', respectively) and the denominator refers to the band or thermally-averaged rate coefficient. Table 6 summarizes the band rate coefficients.

### A.2.3. Collisions with H₂O.

Buffa et al. (2000), Boursier et al. (2020) and Mandal & Babikov (2023) have calculated the rate coefficients for rotational deexcitation within the ground vibrational state of H₂O in collisions with other H₂O molecules. The reported calculations assume that the target molecule sees a Boltzmann distribution of projectiles:

$$H_2O(000, J'_{\tau'}) + H_2O(thermal) \rightarrow H_2O(000, J''_{\tau''}) + H_2O(thermal).$$

This treatment should be adequate here as large deviations of the H₂O states from a Boltzmann distribution only occur at very low densities where the H₂O abundance is likely too small to make a difference in the atmospheric energy budget. In our molecular model, we adopted the Mandal & Babikov (2023) rate coefficients. We complemented them with those by Boursier et al. (2020), as downloaded from BASECOL, which span more states. We used the NM87 method for extrapolation to higher-energy states and for interpolation of a few low-energy states.

The vibrational self-relaxation of H₂O has been investigated experimentally (Kung & Center 1975; Finzi et al. 1977; Zittel & Masturzo 1989, 1991) and theoretically (Shin 1993; Huestis 2006; Braunstein & Conforti 2015). The cross sections and rate coefficients are typically much larger than those for collisions with other heavy particles. There is some evidence that H₂O vibrational self-relaxation occurs through the scheme:

$$ \begin{align} H_2O(010) + X &\rightarrow H_2O(000) + X \nonumber \\ H_2O(020) + X &\rightarrow H_2O(010) + X \tag{A1} \\ H_2O(100,001) + X &\rightarrow H_2O(020) + X \nonumber \end{align} $$

where X≡H₂O. We implemented the above channels using the band rate coefficients reported by Zittel & Masturzo (1989) and omitting their weak temperature dependences. The scheme is consistent with the idea that the two stretching modes, being nearly equal in energy, behave as a single one. For (001)→(100) and for rotational relaxation within each excited vibrational state, we implemented the band rate coefficient for relaxation within the ground vibrational state to ensure rapid thermalization. For these cases, in which the band rate coefficient is known but not the state-to-state partitioning, we assumed that $t_{ji}^{H_2O} = t_{ji}^{H_2}$. Table 7 summarizes this.

The compilation of state-to-state rate coefficients for H, O, O₂ and He is described in the Supplementary Information.
<table><caption>**Table 4.** Radiative properties for the bands connecting some of the lower-energy vibrational states of H₂O. λ<sub>v'v''</sub> is the wavelength. A<sub>v'v''</sub>=∑<sub>j∈v',i∈v'',E<sub>j</sub>>E<sub>i</sub> g<sub>j</sub> exp(-E<sub>j</sub>/kT)A<sub>ji</sub>/∑<sub>j∈v'</sub> g<sub>j</sub> exp(-E<sub>j</sub>/kT) is the transition probability, calculated assuming that the upper states are thermalized at 400 K. g<sub>v''v'</sub><sup>tr1</sup> and g<sub>v''v'</sub><sup>tr1</sup> are the photoexcitation frequencies induced by a Sun- (T<sub>bb</sub>=5,780 K) and a TRAPPIST 1-like (T<sub>bb</sub>=2,500 K) black-body of the true stellar sizes radiating from 1 AU. The photoexcitation frequencies are calculated with Eq. 3 of Crovisier (1984) and assume optically thin conditions. The last two columns represent the energy rates associated with photoexcitation. The tabulated values were calculated for o-H₂O states; the values for p-H₂O states are similar. For preparing this table, we used all rotational states with J≤20 within each vibrational state.</caption><thead><tr><th>v'</th><th>v''</th><th>1/λ<sub>v'v''</sub></th><th>A<sub>v'v''</sub></th><th>g<sub>v''v'</sub><sup>tr1</sup></th><th>g<sub>v''v'</sub><sup>tr1</sup></th><th>g<sub>v''v'</sub>hc/λ<sub>v'v''</sub></th><th>g<sub>v''v'</sub>hc/λ<sub>v'v''</sub></th></tr><tr><th></th><th></th><th>[cm<sup>-1</sup>]</th><th>[s<sup>-1</sup>]</th><th>[s<sup>-1</sup>]</th><th>[s<sup>-1</sup>]</th><th>[erg s<sup>-1</sup>]</th><th>[erg s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>(000)</td><td>(000)</td><td>—</td><td>1.487(+0)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>(010)</td><td>(010)</td><td>—</td><td>1.575(+0)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>(000)</td><td>1594.9</td><td>2.388(+1)</td><td>2.649(-4)</td><td>1.220(-6)</td><td>8.393(-17)</td><td>3.864(-19)</td></tr><tr><td>(020)</td><td>(020)</td><td>—</td><td>1.671(+0)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>(010)</td><td>1556.8</td><td>4.685(+1)</td><td>5.352(-4)</td><td>2.483(-6)</td><td>1.655(-16)</td><td>7.678(-19)</td></tr><tr><td></td><td>(000)</td><td>3151.7</td><td>7.002(-1)</td><td>3.178(-6)</td><td>1.048(-8)</td><td>1.989(-18)</td><td>6.560(-21)</td></tr><tr><td>(100)</td><td>(100)</td><td>—</td><td>1.320(+0)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>(020)</td><td>505.5</td><td>4.905(-2)</td><td>1.978(-6)</td><td>1.116(-8)</td><td>1.986(-19)</td><td>1.121(-21)</td></tr><tr><td></td><td>(010)</td><td>2062.3</td><td>1.246(+0)</td><td>1.004(-5)</td><td>4.203(-8)</td><td>4.113(-18)</td><td>1.722(-20)</td></tr><tr><td></td><td>(000)</td><td>3657.2</td><td>7.701(+0)</td><td>2.804(-5)</td><td>8.212(-8)</td><td>2.037(-17)</td><td>5.965(-20)</td></tr><tr><td>(001)</td><td>(001)</td><td>—</td><td>1.379(+0)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>(100)</td><td>98.8</td><td>1.273(-1)</td><td>2.765(-5)</td><td>1.672(-7)</td><td>5.424(-19)</td><td>3.280(-21)</td></tr><tr><td></td><td>(020)</td><td>604.2</td><td>5.705(-2)</td><td>1.900(-6)</td><td>1.054(-8)</td><td>2.281(-19)</td><td>1.265(-21)</td></tr><tr><td></td><td>(010)</td><td>2161.0</td><td>1.961(+0)</td><td>1.488(-5)</td><td>6.102(-8)</td><td>6.387(-18)</td><td>2.619(-20)</td></tr><tr><td></td><td>(000)</td><td>3755.9</td><td>9.131(+1)</td><td>3.191(-4)</td><td>9.128(-7)</td><td>2.381(-16)</td><td>6.810(-19)</td></tr><tr><td>(030)</td><td>(030)</td><td>—</td><td>1.848(+0)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>(001)</td><td>910.8</td><td>3.383(-2)</td><td>7.187(-7)</td><td>3.771(-9)</td><td>1.300(-19)</td><td>6.823(-22)</td></tr><tr><td></td><td>(100)</td><td>1009.6</td><td>4.857(-2)</td><td>9.191(-7)</td><td>4.736(-9)</td><td>1.843(-19)</td><td>9.498(-22)</td></tr><tr><td></td><td>(020)</td><td>1515.0</td><td>6.912(+1)</td><td>8.158(-4)</td><td>3.816(-6)</td><td>2.455(-16)</td><td>1.148(-18)</td></tr><tr><td></td><td>(010)</td><td>3071.9</td><td>1.696(+0)</td><td>7.985(-6)</td><td>2.681(-8)</td><td>4.872(-18)</td><td>1.636(-20)</td></tr><tr><td></td><td>(000)</td><td>4666.7</td><td>8.912(-3)</td><td>2.195(-8)</td><td>5.009(-11)</td><td>2.035(-20)</td><td>4.643(-23)</td></tr><tr><td>(110)</td><td>(110)</td><td>—</td><td>1.446(+0)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>(030)</td><td>568.5</td><td>3.117(-2)</td><td>1.109(-6)</td><td>6.188(-9)</td><td>1.252(-19)</td><td>6.988(-22)</td></tr><tr><td></td><td>(001)</td><td>1479.3</td><td>9.379(-1)</td><td>1.139(-5)</td><td>5.366(-8)</td><td>3.347(-18)</td><td>1.577(-20)</td></tr><tr><td></td><td>(100)</td><td>1578.1</td><td>2.255(+1)</td><td>2.534(-4)</td><td>1.171(-6)</td><td>7.944(-17)</td><td>3.670(-19)</td></tr><tr><td></td><td>(020)</td><td>2083.5</td><td>2.068(+0)</td><td>1.645(-5)</td><td>6.858(-8)</td><td>6.809(-18)</td><td>2.838(-20)</td></tr><tr><td></td><td>(010)</td><td>3640.3</td><td>6.062(+0)</td><td>2.222(-5)</td><td>6.535(-8)</td><td>1.607(-17)</td><td>4.726(-20)</td></tr><tr><td></td><td>(000)</td><td>5235.2</td><td>6.202(-1)</td><td>1.251(-6)</td><td>2.463(-9)</td><td>1.301(-18)</td><td>2.561(-21)</td></tr><tr><td>(011)</td><td>(011)</td><td>—</td><td>1.488(+0)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>(110)</td><td>96.1</td><td>8.431(-2)</td><td>1.884(-5)</td><td>1.140(-7)</td><td>3.594(-19)</td><td>2.174(-21)</td></tr><tr><td></td><td>(030)</td><td>664.5</td><td>1.886(-2)</td><td>5.668(-7)</td><td>3.110(-9)</td><td>7.482(-20)</td><td>4.105(-22)</td></tr><tr><td></td><td>(001)</td><td>1575.4</td><td>2.113(+1)</td><td>2.380(-4)</td><td>1.100(-6)</td><td>7.447(-17)</td><td>3.442(-19)</td></tr><tr><td></td><td>(100)</td><td>1674.1</td><td>1.280(+0)</td><td>1.339(-5)</td><td>6.069(-8)</td><td>4.453(-18)</td><td>2.018(-20)</td></tr><tr><td></td><td>(020)</td><td>2179.6</td><td>3.690(+0)</td><td>2.770(-5)</td><td>1.131(-7)</td><td>1.199(-17)</td><td>4.899(-20)</td></tr><tr><td></td><td>(010)</td><td>3736.4</td><td>9.916(+1)</td><td>3.493(-4)</td><td>1.004(-6)</td><td>2.593(-16)</td><td>7.451(-19)</td></tr><tr><td></td><td>(000)</td><td>5331.3</td><td>1.850(+1)</td><td>3.611(-5)</td><td>6.931(-8)</td><td>3.824(-17)</td><td>7.340(-20)</td></tr><tr><td>(040)</td><td>(040)</td><td>—</td><td>2.072(+0)</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td></td><td>(011)</td><td>802.6</td><td>8.669(-3)</td><td>2.120(-7)</td><td>1.134(-9)</td><td>3.379(-20)</td><td>1.809(-22)</td></tr><tr><td></td><td>(110)</td><td>898.6</td><td>5.471(-2)</td><td>1.180(-6)</td><td>6.205(-9)</td><td>2.106(-19)</td><td>1.108(-21)</td></tr><tr><td></td><td>(030)</td><td>1467.1</td><td>9.077(+1)</td><td>1.113(-3)</td><td>5.257(-6)</td><td>3.245(-16)</td><td>1.532(-18)</td></tr><tr><td></td><td>(001)</td><td>2377.9</td><td>1.740(-3)</td><td>1.165(-8)</td><td>4.562(-11)</td><td>5.503(-21)</td><td>2.155(-23)</td></tr><tr><td></td><td>(100)</td><td>2476.7</td><td>3.228(-3)</td><td>2.047(-8)</td><td>7.848(-11)</td><td>1.007(-20)</td><td>3.861(-23)</td></tr><tr><td></td><td>(020)</td><td>2982.2</td><td>2.758(+0)</td><td>1.355(-5)</td><td>4.643(-8)</td><td>8.026(-18)</td><td>2.751(-20)</td></tr><tr><td></td><td>(010)</td><td>4539.0</td><td>2.808(-2)</td><td>7.247(-8)</td><td>1.708(-10)</td><td>6.534(-20)</td><td>1.540(-22)</td></tr><tr><td></td><td>(000)</td><td>6133.8</td><td>1.558(-3)</td><td>2.338(-9)</td><td>3.614(-12)</td><td>2.849(-21)</td><td>4.404(-24)</td></tr></tbody></table>
<table><caption>**Table 5.** $k_{v',v''}^Q$: Band rate coefficients for H₂O deexcitation in collisions with thermal electrons. Only the o-H₂O rate coefficients are quoted; those for p-H₂O are similar. Refs.: [FA08], Faure & Josselin (2008); [AY21], Ayouz et al. (2021).</caption><thead><tr><th>v'</th><th>v''</th><th>200 K</th><th>400 K</th><th>800 K</th><th>1200 K</th><th>1600 K</th><th>Source</th></tr></thead><tbody><tr><td>000</td><td>000</td><td>5.08(-07)</td><td>4.97(-07)</td><td>4.59(-07)</td><td>4.23(-07)</td><td>3.96(-07)</td><td>[FA08]</td></tr><tr><td>010</td><td>000</td><td>2.12(-08)</td><td>1.44(-08)</td><td>9.88(-09)</td><td>7.83(-09)</td><td>6.68(-09)</td><td>[FA08], [AY21]</td></tr><tr><td>010</td><td>010</td><td>5.02(-07)</td><td>4.94(-07)</td><td>4.58(-07)</td><td>4.24(-07)</td><td>3.99(-07)</td><td>[FA08]</td></tr><tr><td>020</td><td>000</td><td>7.57(-10)</td><td>5.38(-10)</td><td>3.77(-10)</td><td>3.09(-10)</td><td>2.67(-10)</td><td>[FA08], [AY21]</td></tr><tr><td>020</td><td>010</td><td>3.93(-08)</td><td>2.72(-08)</td><td>1.87(-08)</td><td>1.48(-08)</td><td>1.26(-08)</td><td>[FA08], [AY21]</td></tr><tr><td>020</td><td>020</td><td>4.94(-07)</td><td>4.89(-07)</td><td>4.57(-07)</td><td>4.27(-07)</td><td>4.03(-07)</td><td>[FA08]</td></tr><tr><td>100</td><td>000</td><td>6.94(-09)</td><td>5.03(-09)</td><td>3.70(-09)</td><td>3.09(-09)</td><td>2.74(-09)</td><td>[FA08], [AY21]</td></tr><tr><td>100</td><td>010</td><td>1.32(-08)</td><td>1.19(-08)</td><td>9.24(-09)</td><td>7.70(-09)</td><td>6.68(-09)</td><td>[FA08]</td></tr><tr><td>100</td><td>020</td><td>1.45(-08)</td><td>1.20(-08)</td><td>8.81(-09)</td><td>7.12(-09)</td><td>6.11(-09)</td><td>[FA08]</td></tr><tr><td>100</td><td>100</td><td>5.09(-07)</td><td>4.98(-07)</td><td>4.63(-07)</td><td>4.32(-07)</td><td>4.08(-07)</td><td>[FA08]</td></tr><tr><td>001</td><td>000</td><td>6.63(-10)</td><td>5.90(-10)</td><td>5.63(-10)</td><td>5.50(-10)</td><td>5.42(-10)</td><td>[FA08], [AY21]</td></tr><tr><td>001</td><td>010</td><td>6.72(-10)</td><td>6.36(-10)</td><td>5.12(-10)</td><td>4.74(-10)</td><td>4.59(-10)</td><td>[FA08]</td></tr><tr><td>001</td><td>020</td><td>1.45(-08)</td><td>1.21(-08)</td><td>8.96(-09)</td><td>7.27(-09)</td><td>6.21(-09)</td><td>[FA08]</td></tr><tr><td>001</td><td>100</td><td>4.47(-08)</td><td>3.38(-08)</td><td>2.52(-08)</td><td>2.14(-08)</td><td>1.91(-08)</td><td>[FA08]</td></tr><tr><td>001</td><td>001</td><td>4.94(-07)</td><td>4.86(-07)</td><td>4.54(-07)</td><td>4.25(-07)</td><td>4.03(-07)</td><td>[FA08]</td></tr></tbody></table>

<table><caption>**Table 6.** $k_{v',v''}^Q$: Band rate coefficients for H₂O deexcitation in collisions with a thermalized 3:1-ortho:para population of H₂. Only the o-H₂O rate coefficients are quoted; those for p-H₂O are similar. Refs.: [FA08], Faure & Josselin (2008); [DA11], Daniel et al. (2011).</caption><thead><tr><th>v'</th><th>v''</th><th>200</th><th>400</th><th>800</th><th>1200</th><th>1600</th><th>Source</th></tr></thead><tbody><tr><td>000</td><td>000</td><td>2.04(-10)</td><td>2.54(-10)</td><td>3.41(-10)</td><td>4.24(-10)</td><td>4.88(-10)</td><td>[DA11], [FA08]</td></tr><tr><td>010</td><td>000</td><td>1.29(-12)</td><td>1.55(-12)</td><td>3.56(-12)</td><td>7.12(-12)</td><td>1.23(-11)</td><td>[FA08]</td></tr><tr><td>010</td><td>010</td><td>1.77(-10)</td><td>2.82(-10)</td><td>4.13(-10)</td><td>4.92(-10)</td><td>5.51(-10)</td><td>[FA08]</td></tr><tr><td>020</td><td>000</td><td>5.94(-13)</td><td>7.06(-13)</td><td>1.65(-12)</td><td>3.35(-12)</td><td>5.84(-12)</td><td>[FA08]</td></tr><tr><td>020</td><td>010</td><td>2.29(-12)</td><td>2.74(-12)</td><td>6.34(-12)</td><td>1.28(-11)</td><td>2.25(-11)</td><td>[FA08]</td></tr><tr><td>020</td><td>020</td><td>1.73(-10)</td><td>2.77(-10)</td><td>4.02(-10)</td><td>4.75(-10)</td><td>5.31(-10)</td><td>[FA08]</td></tr><tr><td>100</td><td>000</td><td>7.99(-14)</td><td>9.51(-14)</td><td>2.24(-13)</td><td>4.52(-13)</td><td>7.97(-13)</td><td>[FA08]</td></tr><tr><td>100</td><td>010</td><td>7.99(-14)</td><td>9.51(-14)</td><td>2.24(-13)</td><td>4.52(-13)</td><td>7.97(-13)</td><td>[FA08]</td></tr><tr><td>100</td><td>020</td><td>6.34(-13)</td><td>7.50(-13)</td><td>1.72(-12)</td><td>3.45(-12)</td><td>5.98(-12)</td><td>[FA08]</td></tr><tr><td>100</td><td>100</td><td>1.81(-10)</td><td>2.84(-10)</td><td>4.04(-10)</td><td>4.76(-10)</td><td>5.32(-10)</td><td>[FA08]</td></tr><tr><td>001</td><td>000</td><td>7.95(-14)</td><td>9.58(-14)</td><td>2.30(-13)</td><td>4.71(-13)</td><td>8.45(-13)</td><td>[FA08]</td></tr><tr><td>001</td><td>010</td><td>7.95(-14)</td><td>9.58(-14)</td><td>2.30(-13)</td><td>4.71(-13)</td><td>8.45(-13)</td><td>[FA08]</td></tr><tr><td>001</td><td>020</td><td>6.35(-13)</td><td>7.64(-13)</td><td>1.78(-12)</td><td>3.58(-12)</td><td>6.14(-12)</td><td>[FA08]</td></tr><tr><td>001</td><td>100</td><td>2.85(-11)</td><td>4.05(-11)</td><td>6.46(-11)</td><td>8.70(-11)</td><td>1.03(-10)</td><td>[FA08]</td></tr><tr><td>001</td><td>001</td><td>1.89(-10)</td><td>2.86(-10)</td><td>3.92(-10)</td><td>4.54(-10)</td><td>4.99(-10)</td><td>[FA08]</td></tr></tbody></table>
<table><caption>**Table 7.** $k_{v'v''}^Q$: Band rate coefficients for H₂O deexcitation in collisions with a thermalized 3:1-ortho:para population of H₂O. Only the o-H₂O rate coefficients are quoted; those for p-H₂O are similar. Refs.: [BO20], Boursier et al. (2020); [MA23], Mandal & Babikov (2023); [ZI89], Zittel & Masturzo (1989).</caption><thead><tr><th>v'</th><th>v''</th><th>200</th><th>400</th><th>800</th><th>Source</th></tr></thead><tbody><tr><td>000</td><td>000</td><td>1.03(-09)<br>1.18(-09)</td><td>1.28(-09)<br>1.42(-09)</td><td>1.54(-09)<br>1.80(-09)</td><td>[BO20,MA23]<br>k<sub>u</sub><sup>Q</sup>(↓)</td></tr><tr><td>010</td><td>000</td><td>5.56(-11)</td><td>5.49(-11)</td><td>5.49(-11)</td><td>[ZI89]</td></tr><tr><td>010</td><td>010</td><td>1.02(-09)</td><td>1.28(-09)</td><td>1.54(-09)</td><td>See text</td></tr><tr><td>020</td><td>010</td><td>1.10(-10)</td><td>1.10(-10)</td><td>1.10(-10)</td><td>[ZI89]</td></tr><tr><td>020</td><td>020</td><td>1.03(-09)</td><td>1.28(-09)</td><td>1.54(-09)</td><td>See text</td></tr><tr><td>100</td><td>020</td><td>2.40(-11)</td><td>2.40(-11)</td><td>2.40(-11)</td><td>[ZI89]</td></tr><tr><td>100</td><td>100</td><td>1.03(-09)</td><td>1.28(-09)</td><td>1.54(-09)</td><td>See text</td></tr><tr><td>001</td><td>020</td><td>2.40(-11)</td><td>2.40(-11)</td><td>2.40(-11)</td><td>[ZI89]</td></tr><tr><td>001</td><td>100</td><td>1.03(-09)</td><td>1.28(-09)</td><td>1.54(-09)</td><td>See text</td></tr><tr><td>001</td><td>001</td><td>1.03(-09)</td><td>1.28(-09)</td><td>1.54(-09)</td><td>See text</td></tr></tbody></table>
REFERENCES

Acuña, L., Deleuil, M., Mousis, O., et al. 2021, A&A, 647, A53, doi: 10.1051/0004-6361/202039885

Asensio Ramos, A., & Elitzur, M. 2018, A&A, 616, A131, doi: 10.1051/0004-6361/201731943

Athay, R. G., & Skumanich, A. 1971, ApJ, 170, 605, doi: 10.1086/151246

Ayouz, M., Faure, A., Tennyson, J., Tudorovskaya, M., & Kokoouline, V. 2021, Atoms, 9, 62, doi: 10.3390/atoms9030062

Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087, doi: 10.1111/j.1365-2966.2006.10184.x

Benneke, B., Wong, I., Piaulet, C., et al. 2019, ApJL, 887, L14, doi: 10.3847/2041-8213/ab59dc

Bernstein, L. S., Elgin, J. B., Pike, C. P., et al. 1996, J. Geophys. Res., 101, 383, doi: 10.1029/95JA02893

Bitsch, B., Raymond, S. N., Buchhave, L. A., et al. 2021, A&A, 649, L5, doi: 10.1051/0004-6361/202140793

Boursier, C., Mandal, B., Babikov, D., & Dubernet, M. L. 2020, MNRAS, 498, 5489, doi: 10.1093/mnras/staa2713

Braunstein, M., & Conforti, P. F. 2013, JChPh, 138, 074303, doi: 10.1063/1.4790589

—. 2015, Journal of Physical Chemistry A, 119, 3311, doi: 10.1021/jp5115527

Buffa, G., Tarrini, O., Scappini, F., & Cecchi-Pestellini, C. 2000, ApJS, 128, 597, doi: 10.1086/313389

Cabrera-González, L. D., Denis-Alpizar, O., Páez-Hernández, D., & Stoecklin, T. 2022, MNRAS, 514, 4426, doi: 10.1093/mnras/stac1643

Carone, L., Keppens, R., Decin, L., & Henning, T. 2018, MNRAS, 473, 4672, doi: 10.1093/mnras/stx2732

Cernicharo, J., & Crovisier, J. 2005, SSRv, 119, 29, doi: 10.1007/s11214-005-8058-x

Coppola, C. M., Lique, F., Mazzia, F., Esposito, F., & Kazandjian, M. V. 2019, MNRAS, 486, 1590, doi: 10.1093/mnras/stz927

Crovisier, J. 1984, A&A, 130, 361

Daniel, F., Dubernet, M. L., & Grosjean, A. 2011, A&A, 536, A76, doi: 10.1051/0004-6361/201118049

Daniel, F., Faure, A., Dagdigian, P. J., et al. 2015, MNRAS, 446, 2312, doi: 10.1093/mnras/stu2287

Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nature Astronomy, 5, 775, doi: 10.1038/s41550-021-01381-5

Dubernet, M. L., Daniel, F., Grosjean, A., & Lin, C. Y. 2009, A&A, 497, 911, doi: 10.1051/0004-6361/200810680

Dubernet, M. L., Daniel, F., Grosjean, A., et al. 2006, A&A, 460, 323, doi: 10.1051/0004-6361:20065804

Dunn, M. G., Skinner, G. T., & Treanor, C. E. 1975, AIAA Journal, 13, 803, doi: 10.2514/3.60441

Elitzur, M., & Asensio Ramos, A. 2006, MNRAS, 365, 779, doi: 10.1111/j.1365-2966.2005.09770.x

Faure, A., Crimier, N., Ceccarelli, C., et al. 2007, A&A, 472, 1029, doi: 10.1051/0004-6361:20077678

Faure, A., & Josselin, E. 2008, A&A, 492, 257, doi: 10.1051/0004-6361:200810717

Faure, A., Wiesenfeld, L., Wernli, M., & Valiron, P. 2006, JChPh, 124, 214310, doi: 10.1063/1.2204032

Feofilov, A. G., & Kutepov, A. A. 2012, Surveys in Geophysics, 33, 1231, doi: 10.1007/s10712-012-9204-0

Feofilov, A. G., Kutepov, A. A., Pesnell, W. D., et al. 2009, Atmospheric Chemistry & Physics, 9, 8139, doi: 10.5194/acp-9-8139-2009

Finzi, J., Hovis, F. E., Panfilov, V. N., Hess, P., & Moore, C. B. 1977, JChPh, 67, 4053, doi: 10.1063/1.435379

Funke, B., López-Puertas, M., García-Comas, M., et al. 2012, JQSRT, 113, 1771, doi: 10.1016/j.jqsrt.2012.05.001

Galewsky, J., Steen-Larsen, H. C., Field, R. D., et al. 2016, Reviews of Geophysics, 54, 809, doi: 10.1002/2015RG000512

García Muñoz, A. 2007, Planet. Space Sci., 55, 1426, doi: 10.1016/j.pss.2007.03.007

—. 2023a, Icarus, 392, 115373, doi: 10.1016/j.icarus.2022.115373

—. 2023b, A&A, 672, A77, doi: 10.1051/0004-6361/202245766

García Muñoz, A., Fossati, L., Youngblood, A., et al. 2021, ApJL, 907, L36, doi: 10.3847/2041-8213/abd9b8

García Muñoz, A., Youngblood, A., Fossati, L., et al. 2020, ApJL, 888, L21, doi: 10.3847/2041-8213/ab61ff

García-Vázquez, R., Faure, A., & Stoecklin, T. 2023, ChemPhysChem

Gersch, A. M., Feaga, L. M., & A'Hearn, M. F. 2018, ApJ, 854, 149, doi: 10.3847/1538-4357/aa9795

Gillet, A., García Muñoz, A., & Strugarek, A. 2023, â, 680, A33, doi: 10.1051/0004-6361/202347066

Goldsmith, P. F., & Langer, W. D. 1978, ApJ, 222, 881, doi: 10.1086/156206

Gray, M. D., Baudry, A., Richards, A. M. S., et al. 2016, MNRAS, 456, 374, doi: 10.1093/mnras/stv2437

Green, S., Maluendes, S., & McLean, A. D. 1993, ApJS, 85, 181, doi: 10.1086/191760

Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555, doi: 10.1086/190631

Hrodmarsson, H. R., & van Dishoeck, E. F. 2023, A&A, 675, A25, doi: 10.1051/0004-6361/202346645

Huestis, D. L. 2006, Journal of Physical Chemistry A, 110, 6638, doi: 10.1021/jp054889n
Ih, J., Kempton, E. M. R., Whittaker, E. A., & Lessard, M. 2023, ApJL, 952, L4, doi: 10.3847/2041-8213/ace03b

Irons, F. E. 1978, MNRAS, 182, 705, doi: 10.1093/mnras/182.4.705

Izidoro, A., & Piani, L. 2023, arXiv e-prints, arXiv:2302.02674, doi: 10.48550/arXiv.2302.02674

Johnson, B. R. 1986, JChPh, 84, 176, doi: 10.1063/1.450807

Kaspi, Y., & Showman, A. P. 2015, ApJ, 804, 60, doi: 10.1088/0004-637X/804/1/60

Kimura, T., & Ikoma, M. 2022, Nature Astronomy, 6, 1296, doi: 10.1038/s41550-022-01781-1

Kite, E. S., & Schaefer, L. 2021, ApJL, 909, L22, doi: 10.3847/2041-8213/abe7dc

Kolb, C. E., & Elgin, J. B. 1977, JChPh, 66, 119, doi: 10.1063/1.433658

Kung, R. T. V., & Center, R. E. 1975, JChPh, 62, 2187, doi: 10.1063/1.430786

Lang, B., Breitegger, P., Brunnhofer, G., et al. 2020, Applied Physics B: Lasers and Optics, 126, 64, doi: 10.1007/s00340-020-7409-3

Lincowski, A. P., Meadows, V. S., Zieba, S., et al. 2023, ApJL, 955, L7, doi: 10.3847/2041-8213/acee02

López-Morales, M., Ben-Ami, S., Gonzalez-Abad, G., et al. 2019, AJ, 158, 24, doi: 10.3847/1538-3881/ab21d7

López-Puertas, M., & Taylor, F. W. 2001, Non-LTE radiative transfer in the atmosphere, Vol. 3

Madhusudhan, N., Nixon, M. C., Welbanks, L., Piette, A. A., & Booth, R. A. 2020, ApJL, 891, L7, doi: 10.3847/2041-8213/ab7229

Mandal, B., & Babikov, D. 2023, A&A, 671, A51, doi: 10.1051/0004-6361/202245699

Manuilova, R. O., Feofilov, A. G., Kutepov, A. A., & Yankovsky, V. A. 2015, Advances in Space Research, 56, 1806, doi: 10.1016/j.asr.2014.12.002

Meyerott, R. E., Swenson, G. R., Schweitzer, E. L., & Koch, D. G. 1994, J. Geophys. Res., 99, 17559, doi: 10.1029/94JA01241

Moran, S. E., Stevenson, K. B., Sing, D. K., et al. 2023, ApJL, 948, L11, doi: 10.3847/2041-8213/accb9c

Morris, M. A., Desch, S. J., & Ciesla, F. J. 2009, ApJ, 691, 320, doi: 10.1088/0004-637X/691/1/320

Neufeld, D. A., & Kaufman, M. J. 1993, ApJ, 418, 263, doi: 10.1086/173388

Neufeld, D. A., & Melnick, G. J. 1987, ApJ, 322, 266, doi: 10.1086/165721

Phillips, T. R., Maluendes, S., & Green, S. 1996, ApJS, 107, 467, doi: 10.1086/192372

Piaulet, C., Benneke, B., Almenara, J. M., et al. 2023, Nature Astronomy, 7, 206, doi: 10.1038/s41550-022-01835-4

Redmon, M. J., Schatz, G. C., & Garrett, B. C. 1986, JChPh, 84, 764, doi: 10.1063/1.450574

Sánchez-Lavega, A., Irwin, P., & García Muñoz, A. 2023, A&A Rv, 31, 5, doi: 10.1007/s00159-023-00150-9

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369, doi: 10.1051/0004-6361:20041729

Shin, H. K. 1993, JChPh, 98, 1964, doi: 10.1063/1.464230

Smith, P. L., Yoshino, K., Griesinger, H. E., & Black, J. H. 1981, ApJ, 250, 166, doi: 10.1086/159359

Snellen, I. 2014, Philosophical Transactions of the Royal Society of London Series A, 372, 20130075, doi: 10.1098/rsta.2013.0075

Stoecklin, T., Cabrera-González, L. D., Denis-Alpizar, O., & Páez-Hernández, D. 2021, JChPh, 154, 144307, doi: 10.1063/5.0047718

Stoecklin, T., Denis-Alpizar, O., Clergerie, A., et al. 2019, Journal of Physical Chemistry A, 123, 5704, doi: 10.1021/acs.jpca.9b04052

Traub, W. A., & Oppenheimer, B. R. 2010, in Exoplanets, ed. S. Seager, 111–156

Valencia, D., Sasselov, D. D., & O’Connell, R. J. 2007, ApJ, 656, 545, doi: 10.1086/509800

van Dishoeck, E. F., Herbst, E., & Neufeld, D. A. 2013, Chemical Reviews, 113, 9043, doi: 10.1021/cr4003177

Venturini, J., Guilera, O. M., Haldemann, J., Ronco, M. P., & Mordasini, C. 2020, A&A, 643, L1, doi: 10.1051/0004-6361/202039141

Villanueva, G. L., Hammel, H. B., Milam, S. N., et al. 2023, Nature Astronomy, doi: 10.1038/s41550-023-02009-6

Walker, K. M., Yang, B. H., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2014, ApJ, 790, 96, doi: 10.1088/0004-637X/790/2/96

Wiesenfeld, L. 2021, JChPh, 155, 071104, doi: 10.1063/5.0058755

—. 2022, JChPh, 157, 174304, doi: 10.1063/5.0102279

Wordsworth, R. 2015, ApJ, 806, 180, doi: 10.1088/0004-637X/806/2/180

Xie, X., & Mumma, M. J. 1992, ApJ, 386, 720, doi: 10.1086/171053

Yang, B., Nagao, M., Satomi, W., Kimura, M., & Stancil, P. C. 2013, ApJ, 765, 77, doi: 10.1088/0004-637X/765/2/77

Yankovsky, V., Manuilova, R., Babaev, A., Feofilov, A., & Kutepov, A. 2011, International Journal of Remote Sensing, 32, 3065, doi: 10.1080/01431161.2010.541506

Yun, Y. J., Park, Y. S., & Lee, S. H. 2009, A&A, 507, 1785, doi: 10.1051/0004-6361/200912800
Zel'dovich, Y. B., & Raizer, Y. P. 2002, Physics of shock

waves and high-temperature hydrodynamic phenomena

(Dover Publications, Inc.)

Zhou, D. K., Pendleton, W. R., J., Bingham, G. E., Steed, A. J., & Dean, D. A. 1994, Geophys. Res. Lett., 21, 613, doi: 10.1029/94GL00199

Zittel, P. F., & Masturzo, D. E. 1989, JChPh, 90, 977, doi: 10.1063/1.456122

—. 1991, JChPh, 95, 8005, doi: 10.1063/1.461331

Żóltowski, M., Lique, F., Karska, A., & Żuchowski, P. S. 2021, MNRAS, 502, 5356, doi: 10.1093/mnras/stab453
# SUPPLEMENTARY INFORMATION. I. SOME NOTES ON MOLPOP-CEP

For our work, we have downloaded and used version 51eef94 of MOLPOP-CEP (Asensio Ramos & Elitzur 2018) available on the site in Jan 2023. Our approach has been to make the minimum of changes needed. They are of two general types:

* We have updated the molecular model (states, collisional and radiative properties) of H₂O as described elsewhere in this document. We are making publicly available the corresponding data files.

* We have modified 3 routines of MOLPOP-CEP to produce the outputs needed in the calculation of the net cooling rate and other ancillary information useful to gain insight into the solutions. The specific modifications are described in detail below for each of the routines. We are making publicly available the modified MOLPOP-CEP files, the input files to run the code and reproduce our solutions and the output files. The new outputs are written out into the four files: 'steam.lines.dat'; 'steam.coolingrates.dat'; 'steam.populations.dat'; 'steam.transitions.dat'.

Our simulations use the Coupled Escape Probability formalism, which produces exact solutions to the NLTE problem in a static medium. We use the code option to deal with media of spatially-varying properties and omit the effects of dust. We have optionally considered the effects of external irradiation in the simulations, which are implemented through boundary conditions on the 'left' (= bottom of the atmosphere) and 'right' (= top of the atmosphere) ends of the model. We have set a convergence threshold defined by 'Accuracy in solution of the equations = 1.0e-4', which means that the maximum relative difference between two consecutive iterations for any molecular state and altitude in the atmosphere is less than the prescribed threshold. This gives undue weight to the negligibly populated states, which are often responsible for the overall convergence rate, but we accept this treatment because it provides extra confidence in the convergence of the solutions. Additional evidence that the solutions have converged comes from the fact that the cooling rates obtained from Eqs. 7 and 10 generally match to within 3-4 significant digits, which would be difficult to explain if the numerical solution was not fully converged. Lastly, we ran a few additional simulations at all three kinetic temperatures of 200, 800 and 1,500 K with a convergence threshold of '1.0e-8'. The results obtained with both convergence thresholds matched to 3-4 significant digits. We see no reasons to question the convergence of the numerical solutions presented here.

## CEP_MOLPOP_INTERFACE.F90

The routine produces one new output file:

* steam.transitions.dat. It reports, for each pair of upper and lower states and location in the spatial grid, the transition probability and the deexcitation collisional rate of the background gas.

## ESCAPE_CEP.F90

The routine produces one new output file:

* steam.lines.dat. It reports, for each radiative transition and location in the spatial grid, information about the lines such as: states, wavelength, number density of the upper states, net radiative bracket, line emission.

## IO_CEP.F90

The routine produces two new output files:

* steam.populations.dat. It reports the number densities of all the H₂O states at each location in the spatial grid. Both NLTE and LTE densities are reported.

* steam.coolingrates.dat. It reports the cooling rates at each location in the spatial grid. Various forms of the cooling rate are reported, including the net cooling rate $\Gamma^{\text{BB}}$ as calculated through both Eqs. 7 and 10 and the vibrational components $\Gamma_v^{\text{BB}}$.

Upon exit, the routine writes out the transition probabilities $A_{v'v''}$ and rate coefficients for deexcitation between vibrational states $k_{v'v''}^Q$. They are calculated at each location in the spatial grid by adding the contributions from the states within the upper vibrational states $v'$ as if they were in LTE.
# SUPPLEMENTARY INFORMATION. II. COLLISIONS WITH OTHER HEAVY PARTICLES

## COLLISIONS WITH H.

Daniel et al. (2015) have calculated the rate coefficients for rotational deexcitation within the H₂O ground vibrational state in collisions with the H atom. We adopted them from BASECOL and extrapolated them to higher-energy states with the NM87 method. For rotational deexcitation within the excited vibrational state (010), we used the branching ratios $t_{ji}^{H_2}$ and assumed that the band rate coefficient is the same as for rotational relaxation within (000). For vibrational deexcitation (010)→(000), we used the branching ratios $t_{ji}^{H_2}$ and estimated the band rate coefficient from Fig. 9 (bottom; temperatures above 100 K) of Cabrera-González et al. (2022). Table 1 summarizes this.

## COLLISIONS WITH O.

To our knowledge, the only determinations of the H₂O deexcitation rate coefficients in collisions with O atoms are the measurements for states (100,001) by Zittel & Masturzo (1989), typically assumed to relax into (020) (Funke et al. 2012). There have been various investigations of the cross sections for rotational excitation within the ground vibrational state and for excitation from this into other vibrational states (Dunn et al. 1975; Kolb & Elgin 1977; Johnson 1986; Redmon et al. 1986; Meyerott et al. 1994; Zhou et al. 1994; Bernstein et al. 1996; Braunstein & Conforti 2013). These investigations have focused on high collision energies, and it is not clear how to derive from them rate coefficients at temperatures <1,600 K. To fill the gap for the purely rotational channels:

$$H_2O(000) + O \rightarrow H_2O(000) + O$$

$$H_2O(010) + O \rightarrow H_2O(010) + O$$

and the vibrational channels:

$$H_2O(010) + O \rightarrow H_2O(000) + O$$

$$H_2O(020) + O \rightarrow H_2O(010) + O$$

we adopted the corresponding rate coefficients for collisions with H₂ and scaled them by $(\mu_{H_2O-H_2}/\mu_{H_2O-O})^{1/2} \approx 0.46$. Table 2 summarizes this.

## COLLISIONS WITH O₂.

To our knowledge, there are no determinations of the rate coefficients for H₂O rotational deexcitation within the ground vibrational state in collisions with the O₂ molecule. We filled this gap by adopting the rate coefficients for collisions with H₂ and scaling them by $(\mu_{H_2O-H_2}/\mu_{H_2O-O_2})^{1/2} \approx 0.40$. We proceeded similarly for the rotational relaxation within the vibrational state (010).

The collisions with O₂ that involve a change in the H₂O vibrational state are of interest in the NLTE modelling of air (Feofilov et al. 2009; Yankovsky et al. 2011; Funke et al. 2012; Manuilova et al. 2015; Lang et al. 2020). The usual scheme is similar to that of Eq. A1 with $X \equiv O_2(0)$. The relaxation rate coefficient for states (100,001) has been determined by Finzi et al. (1977). The rate coefficient for H₂O(010)+O₂(0)→H₂O(000)+O₂(0) is unclear but is expected to be much lower than that for the nearly-resonant channel H₂O(010)+O₂(0)↔H₂O(000)+O₂(1) (Huestis 2006). Similarly, H₂O(020)+O₂(0)→H₂O(010)+O₂(0) is expected to be slower than H₂O(020)+O₂(0)↔H₂O(010)+O₂(1). This entails that the H₂O and O₂ NLTE problems are coupled and in principle must be treated together, as is often done in the investigation of the Earth's atmosphere. We compiled the rate coefficients listed in Table 3, from which it is possible to assess the significance of collisions with O₂ in the H₂O NLTE problem. For example, Tables 6-7 show that the relaxation of H₂O(010) in collisions with H₂ and H₂O proceeds with rate coefficients on the order of a few×10⁻¹² and 5.5×10⁻¹¹ cm³s⁻¹, respectively, whereas the rate coefficient for collisions with O₂(0) (proceeding through the nearly-resonant channel) is ~7×10⁻¹³ cm³s⁻¹ at 500 K (Huestis 2006, fig. 3). The latter channel appears to be notably slower than the others and suggests that unless O₂ is much more abundant than the other particles, the contribution of O₂ to the H₂O NLTE problem should be minor.
## COLLISIONS WITH HE.

Green et al. (1993) and Yang et al. (2013) have calculated the rate coefficients for H₂O rotational deexcitation within the ground vibrational state in collisions with He atoms. We adopted the less comprehensive but more recent rate coefficients of Yang et al. (2013) and complemented them with those of Green et al. (1993). We extrapolated the resulting set to higher-energy states with the NM87 method. To our knowledge, there are no calculations of rotational deexcitation rates within the vibrational state (010), and we estimated them by scaling the rate coefficients for collisions with H₂ by $(\mu_{H_2O-H_2}/\mu_{H_2O-He})^{1/2} \approx 0.74$.

Stoecklin et al. (2021) have calculated the (010)→(000) rate coefficients, which compare acceptably well with the high-temperature measurements of Kung & Center (1975). We combined the thermally-averaged rate coefficients of Stoecklin et al. (2021) (their fig. 7, bottom) with the branching ratios $t_{ji}^{H_2}$ to estimate the state-to-state (010)→(000) rate coefficients. Table 4 summarizes this.

<table><caption>**Table 1.** $k_{v'v''}^Q$: Band rate coefficients for deexcitation of H₂O in collisions with H atoms. Only the o-H₂O rate coefficients are quoted; those for p-H₂O are similar. Refs.: [DA15], Daniel et al. (2015); [CA22], Cabrera-González et al. (2022).</caption><thead><tr><th><i>v'</i></th><th><i>v''</i></th><th>200</th><th>400</th><th>800</th><th>1200</th><th>1500</th><th>Source</th></tr></thead><tbody><tr><td>000</td><td>000</td><td>3.95(-11)</td><td>8.38(-11)</td><td>1.61(-10)</td><td>2.21(-10)</td><td>2.59(-10)</td><td>[DA15]</td></tr><tr><td></td><td></td><td>7.10(-11)</td><td>1.14(-10)</td><td>1.89(-10)</td><td>2.46(-10)</td><td>2.82(-10)</td><td>$k_{ii}^Q(\downarrow)$</td></tr><tr><td>010</td><td>000</td><td>5.93(-14)</td><td>1.76(-13)</td><td>5.24(-13)</td><td>9.89(-13)</td><td>1.56(-12)</td><td>[CA22]</td></tr><tr><td>010</td><td>010</td><td>3.96(-11)</td><td>8.37(-11)</td><td>1.61(-10)</td><td>2.21(-10)</td><td>2.57(-10)</td><td>See text</td></tr></tbody></table>

<table><caption>**Table 2.** $k_{v'v''}^Q$: Band rate coefficients for H₂O deexcitation in collisions with O atoms. Refs.: [ZI89], Zittel & Masturzo (1989); [FU12], Funke et al. (2012).</caption><thead><tr><th><i>v'</i></th><th><i>v''</i></th><th>200</th><th>400</th><th>800</th><th>1200</th><th>1600</th><th>Source</th></tr></thead><tbody><tr><td>000</td><td>000</td><td>9.38(-11)</td><td>1.17(-10)</td><td>1.57(-10)</td><td>1.95(-10)</td><td>2.24(-10)</td><td>μ scaling</td></tr><tr><td>010</td><td>000</td><td>5.93(-13)</td><td>7.13(-13)</td><td>1.64(-12)</td><td>3.27(-12)</td><td>5.66(-12)</td><td>μ scaling</td></tr><tr><td>010</td><td>010</td><td>8.14(-11)</td><td>1.30(-10)</td><td>1.90(-10)</td><td>2.26(-10)</td><td>2.53(-10)</td><td>μ scaling</td></tr><tr><td>020</td><td>010</td><td>1.05(-12)</td><td>1.26(-12)</td><td>2.92(-12)</td><td>5.89(-12)</td><td>1.03(-11)</td><td>μ scaling</td></tr><tr><td>100,001</td><td>020</td><td>2.45(-12)</td><td>3.46(-12)</td><td>4.90(-12)</td><td>6.00(-12)</td><td>6.93(-12)</td><td>[ZI89,FU12]</td></tr></tbody></table>

<table><caption>**Table 3.** $k_{v'v''}^Q$: Band rate coefficients for deexcitation of H₂O in collisions with O₂. References: [FI77], Finzi et al. (1977); [HU06], Huestis (2006); [FU12], Funke et al. (2012). †: Rate coefficients actually refer to the nearly-resonant channels H₂O(010)+O₂(0)→H₂O(000)+O₂(1) and H₂O(020)+O₂(0)→H₂O(010)+O₂(1).</caption><thead><tr><th><i>v'</i></th><th><i>v''</i></th><th>200</th><th>400</th><th>800</th><th>1200</th><th>1600</th><th>Source</th></tr></thead><tbody><tr><td>000</td><td>000</td><td>8.16(-11)</td><td>1.02(-10)</td><td>1.36(-10)</td><td>1.69(-10)</td><td>1.95(-10)</td><td>μ scaling</td></tr><tr><td>010</td><td>000</td><td>7.00(-13)</td><td>7.00(-13)</td><td>7.00(-13)</td><td>7.00(-13)</td><td>7.00(-13)</td><td>[HU06]†</td></tr><tr><td>010</td><td>010</td><td>7.08(-11)</td><td>1.13(-10)</td><td>1.65(-10)</td><td>1.97(-10)</td><td>2.20(-10)</td><td>μ scaling</td></tr><tr><td>020</td><td>010</td><td>2.00(-12)</td><td>2.00(-12)</td><td>2.00(-12)</td><td>2.00(-12)</td><td>2.00(-12)</td><td>[FU12]†</td></tr><tr><td>100,001</td><td>020</td><td>3.30(-13)</td><td>3.30(-13)</td><td>3.30(-13)</td><td>3.30(-13)</td><td>3.30(-13)</td><td>[FI77]</td></tr></tbody></table>
<table><caption>**Table 4.** $k_{v',v''}^Q$: Band rate coefficients for H₂O deexcitation in collisions with He atoms. Only the o-H₂O rate coefficients are reported; those for p-H₂O are similar. Refs.: [GR93], Green et al. (1993); YA13, Yang et al. (2013); ST21, Stoecklin et al. (2021).</caption><thead><tr><th>v'</th><th>v''</th><th>200 K</th><th>400 K</th><th>800 K</th><th>Source</th></tr></thead><tbody><tr><td>000</td><td>000</td><td>4.21(-11)</td><td>7.48(-11)</td><td>1.27(-10)</td><td>[GR93], [YA13]</td></tr><tr><td></td><td></td><td>3.51(-11)</td><td>7.24(-11)</td><td>1.37(-10)</td><td>$k_{u''}^Q(\downarrow)$</td></tr><tr><td>010</td><td>000</td><td>1.81(-15)</td><td>1.38(-14)</td><td>1.00(-13)</td><td>[ST21]</td></tr><tr><td>010</td><td>010</td><td>1.31(-10)</td><td>2.09(-10)</td><td>3.06(-10)</td><td>$\mu$ scaling</td></tr></tbody></table>