© The Authors 2024

1 Introduction

Water is ubiquitous in the Universe, both in the gas and solid ice phases. In particular, the detection of water vapour in astronomical environments, with high kinetic and radiation temperatures, has led to the computation of extensive line lists for ro-vibrational transitions. The most recent one was published in Polyansky et al. (2018) and is available from the ExoMol database¹. This line list is based on available experimental data, comprising all energy levels up to 41 000 cm⁻¹ and a rotational angular momentum up to j = 72. It can be used to model exoplanetary and cometary atmospheres, star-forming regions, and also hot objects such as sunspots, dwarf stars, and giant stars, where temperatures can exceed 3000 K (see Polyansky et al. (2018) and references therein). In these environments, deviations from local thermodynamic equilibrium (LTE) can be substantial, as exemplified by H₂O maser emissions from vibrationally excited states in evolved stars Gray et al. (2016); Baudry et al. (2023). In addition to radiative rates, the collisional rate coefficients are thus critical parameters required for accurately modeling the ro-vibrational excitation of H₂O in space environments.

In the interstellar (ISM) and circumstellar (CSM) media, collisional excitation may be due to atomic and molecular hydrogen, helium atoms, and free electrons. The rotational excitation of H₂O in its ground vibrational state ((000) in normal mode notation) by neutrals has been widely studied, both theoretically and experimentally (see Daniel et al. 2015; Bergeat et al. 2022 and references therein). In contrast, ro-vibrational calculations are scarce and restricted to the (de)excitation of the first excited bending mode (010) (see Wiesenfeld 2021; Garcia-Vazquez et al. 2024 and references therein). In the case of water-electron collisions, both rotational and vibrational excitations have been largely studied, as reviewed in Song et al. (2021). Until recently, however, calculations and measurements were limited to the excitation of the very low-lying vibrational states, namely (010), (100) and (001), corresponding to one-quantum excitations.

In 2021, Ayouz et al. (2021) have considered for the first time both one-quantum and two-quanta excitations in the vibrational excitation of H₂O by electrons. This work was based on a theoretical approach combining the normal mode approximation, a vibrational frame transformation, and the UK R-matrix codes Tennyson (2010); Tennyson et al. (2007); Carr et al. (2012). The overall agreement with measured cross sections was found to be within experimental uncertainties for the bending mode, (000) → (010), and within about a factor of 2 for the stretching modes, (000) → (100) + (001). cross sections for these one-quantum transitions were found to be larger than those for two-quanta transitions by factors of 3–30. Thermally averaged rate coefficients in the range 10–10000 K were also derived from the cross sections and they were fitted to simple analytical formula for use in models.

In the present work, we further extend the calculations by Ayouz et al. (2021) to treat all normal modes simultaneously so that cross sections for inter-mode transitions can be also evaluated. Moreover, the 13 lowest vibrational levels of H₂O are considered, namely, up to level (101), which lies ~7200 cm⁻¹ above the ground state level. The collisional data for such high-lying vibrational states are critically needed to model hot astrophysical environments. For example, in oxygen-rich late-type giant and supergiant stars, rotational transitions in vibrational states up to (011) have been recently detected with the Atacama Large Millimeter/submillimeter Array (ALMA) Baudry et al. (2023).

This article is organized as follows. The next section presents the theoretical approach used in the calculations. In Sect. 3, the obtained rate coefficients for vibrational (de-)excitation are discussed and compared with the data available in literature. In Sect. 4, we assess the uncertainties of the obtained theoretical results. Section 5 presents the conclusions of this study.

2 Theoretical approach

Water is a closed-shell molecule, having the symmetry of the C_2υ point group at equilibrium and the ground state electronic configuration of:

$${\text{X}}^1{A}_1:1a_1^{2}2a_1^{2}1b_2^{2}3a_1^{2}1b_1^2.$$

Water (H₂O) is characterized by three normal modes of vibration: bending, symmetric stretching, and asymmetric stretching with respective frequencies of ω₂, ω₁ and ω₃. In the discussion below, we use dimensionless normal coordinates q₂, q₁, and q₃, corresponding to the three modes.

2.1 Fixed-geometry scattering matrix

The frequencies of the normal modes and the matrix of transformation between the normal modes and Cartesian coordinates required for the present calculations were determined using the ab initio quantum chemistry package MOLPRO Werner et al. (2012). Details on that calculations are given in Table 1 of our previous study by Ayouz et al. (2021).

As a first step in the calculations, the reactance matrix K(q) (K-matrix) is obtained numerically for seven fixed geometry configurations, q, of the molecule, listed in Table 1. For these calculations, the UK R-matrix codes Tennyson (2010); Carr et al. (2012) with the Quantemol-N interface Tennyson et al. (2007) were employed. Channels (indexes) of the matrix are numerated by the angular momentum, l, of the electron and its projection, λ, on the axis perpendicular to the plane of the molecule.

For each molecular geometry, q, in the table above, K(q) is obtained in the C₁ point group with the molecule being in its ground electronic state. The two 1a² core electrons are frozen and eight electrons are kept distributed in the active space including 2 – 9a molecular orbitals. The R-matrix sphere of radius of 10 bohrs and a partial-wave expansion with continuum Gaussian-type orbitals up to l ≤ 4 were used. This calculation with the cc-pVTZ basis set and the described complete active space (CAS₁) are referred to as model 1 hereafter. The other models used to assess the uncertainty of the obtained results are introduced below.

The K-matrices were transformed into scattering matrix S (q) and then the vibrational frame transformation was performed to obtain the scattering matrix in the representation of vibrational states of the target molecule.

Fig. 1 Employed normal coordinate grid (top panel) to evaluate the integral in Eq. (3). In the lower panel, the corresponding bond coordinates grid is represented.

2.2 Vibrational dynamics

Constructing elements of the scattering matrix for transitions from one vibrational state (υ₁, υ₂, υ₃) = υ to another (υ′₁, υ′₂, υ′₃) = υ′, one needs vibrational wave functions. For this purpose the Schrödinger equation for vibrational motion along q,

$$\left[{\displaystyle \sum _{i=1}^3-}\frac{\hslash {\omega }_i}{2}\frac{{\partial }^2}{\partial q_i^2}+V(q)\right]{\psi }_v(q)={ϵ}_v{\psi }_v(q)$$(1)

is solved numerically using a DVR-type method Kokoouline et al. (1999). In the equation, V(q) is the potential energy surface (PES) of H₂O. In the numerical calculations, a grid of 10 × 10 × 10 points was used with q₁, q₂ and q₃ varying from −4 to 4 (dimensionless). For these normal-mode coordinate intervals, intervals for bound coordinates r₁, r₂, θ of H₂O are from 0.47 Å to 1.55 Å for r₁ and r₂, and from 51.92° to 169.89° for θ. Figure 1 displays both grids.

In Eq. (1), we employed the H₂O PES of Mizus et al. (2018) obtained in an ab initio calculations and refined using empirical rovibrational energy levels. This PES is given in terms of bond coordinates (r₁, r₂, θ). It is converted to the normal coordinates (q) using the transformation matrix and the frequencies. The lower panel of Fig. 2 shows the minimum energy path (MEP) to dissociation along r₁ : for each value of r₁, the minimum with respect to r₂ and θ is found and plotted in the figure. The top panel gives the values of r₂ and θ corresponding to the MEP. The global minimum of the PES, 0.1002 eV, is situated at r₁ = r₂ = 0.9578 Å and θ = 104.5101°. Table 2 gives energies of the lowest vibrational states obtained with the PES solving Eq. (1) and compares the energies with the available experimental data.

Fig. 2 Minimum energy path of H2O (lower panel) in eV as a function of r1 and the second bond length r2 and the bond angle θ (top panel). The bond length, r2, is given with the solid line and with values indicated on the left axis, while the bond angle is given with the dashed line with values on the right axis of the top panel. The vertical dashed line represents the position of the global minimum.

2.3 Scattering matrix accounting for vibrational motion

Having vibrational wave functions and geometry-fixed S -matrix depending on the normal coordinates, the scattering matrix describing the process:

$$e^{-}(l\lambda )+{\text{H}}_2\text{O}(v)\to e^{-}\left(l^{\prime }{\lambda }^{\prime }\right)+{\text{H}}_2\text{O}\left(v^{\prime }\right),$$(2)

expressed as

$${\mathcal{𝒮}}_{l^{\prime }{\lambda }^{\prime }{\nu }^{\prime },l\lambda v}^{\text{ex}}=\langle {\psi }_{v^{\prime }}(q)\left|S_{l^{\prime }{\lambda }^{\prime },l\lambda }^{\text{ex}}(q)\right|{\psi }_v(q)\rangle ,$$(3)

where the brackets imply an integration over the three normal coordinates. The integral is evaluated on DVR grid points. The geometry-fixed scattering matrix $S_{l^{\prime }{\lambda }^{\prime },l\lambda }^{\text{ex}}$(q) is obtained in the following way for each scattering energy E_el. First, K(q) is calculated for geometries given in Table 1. Then, the Eigenvalues k(q) of K(q) and the matrix U of eigenvectors are computed. The eigenvalues were used to obtain eigen phaseshifts δ(q) = arctan (k(q)). Next, the matrix of scattering phaseshifts Δ(q) is then built as Δ(q) = U δ(q) U⁻¹, where δ(q) a diagonal matrix made of the eigen phaseshifts δ(q). Finally, a fitting procedure was performed for each matrix element ${\text{Δ}}_{l^{\prime }{\lambda }^{\prime },l\lambda }(q)$ using an analytical expression. Depending on the partial wave indexes, there are two types of elements that are transforming differently, with respect to the change of sign of the asymmetric stretching coordinate q₃: the elements that change sign after the transformation (antisymmetric with respect to the q₃ → −q₃ transformation) and the elements that don't change sign (symmetric under q₃ → −q₃). For the antisymmetric elements, the fitting formula includes only odd powers of q₃ (linear and cubic):

$${\text{Δ}}_{{\mathcal{A}}^{\prime }{l}^{\prime },\lambda l}(q)=c_{0,0}+c_{1,1}{q}_1+c_{2,1}{q}_2+c_{3,1}{q}_3+c_{3,3}{q}_3^3,$$(4)

while for the symmetric elements, the formula includes only even powers of q₃ ($q_3^2$ and $q_3^4$):

$${\text{Δ}}_{{\lambda }^{\prime }{l}^{\prime },\lambda l}(q)=c_{0,0}+c_{1,1}{q}_1+c_{2,1}{q}_2+c_{3,2}{q}_3^2+c_{3,4}{q}_3^4.$$(5)

The fitting coefficients c_{i, j} are obtained using the values of the matrix elements computed for the geometries listed in Table 1. The matrix of phaseshifts is fitted to the analytical forms of Eqs. (4) and (5) for all scattering energies E_el, varying from 10 meV to 6 eV with a step of 10 meV, for which the R-matrix calculations were performed. The obtained fitted form of the matrix will be referred as ${\text{Δ}}_{l^{\prime },{\lambda }^{\prime },l\lambda }^{\text{ex}}(q)$ to distinguish it from the matrix of Eqs. (4) and (5), obtained at the grid of geometries. The scattering matrix, S^ex(q), is deduced from Δ^ex(q) in a similar manner: diagonalizing Δ^ex(q) and constructing S^ex(q) = U^ex exp [2iδ^ex(q)] (U^ex)⁻¹ (for each electron scattering energy), where δ^ex(q) is the diagonal matrix of the eigen phaseshifts.

The elements of the matrix ${\mathcal{S}}_{l^{\prime }{\lambda }^{\prime }{v}^{\prime },l\lambda v}^{\text{e}}$ represent scattering amplitudes for transition from one vibrational-electronic level (lλv) to another (l′λ′v′). In its fitted form, the scattering matrix S^ex(q) can be evaluated at any geometry. It is needed to perform the vibrational frame transformation of Eq. (3) numerically. The vibrational wavefuctions are depicted in Figs. 3–4.

Fig. 3 Wave functions of several lowest vibrational levels (below 7000 cm−1) of H2O. Wave functions are represented as intensity plots depending on the two (out of the three) coordinates while the third one fixed at the value 0.8. On the top of each panel, the vibrational levels and normal coordinates are indicated.

3 Results

3.1 Cross sections

The cros section for the v′ ← v process is obtained from the above scattering matrix

$${\sigma }_{{\nu }^{\prime }\leftarrow v}\left(E_{\text{el}}\right)=\frac{\pi {\hslash }^2}{2m{E}_{\text{el}}}{\displaystyle \sum _{l^{\prime }{\lambda }^{\prime }l\lambda }{\left|S_{l^{\prime }{\lambda }^{\prime }{v}^{\prime },l\lambda v}^{\text{ex}}-{\delta }_{v^{\prime },v}{\delta }_{l^{\prime }{\lambda }^{\prime },l\lambda }\right|}^2},$$(6)

where m and E_el are the mass and energy of the incident electron.

Figures 5 and 6 present a sample of cross sections for the initial states (000) and (110). We can observe that the largest cross sections are those for one-quantum (de)excitations but that inter-mode transitions are also significant. Comparisons to available experimental data and previous calculations, as in Ayouz et al. (2021), can be found in Fig. 7. The cross section data have been added as supplementary material.

Fig. 4 Figure details similar to Fig. 3, with several other wave functions.

Fig. 5 Cross sections for vibrational excitation from the vibrational level (000) to several other levels (υ′1υ′2υ′3) shown on the left. H2O low-lying vibrational levels are depicted with the same colour code on the right.

Fig. 6 Cross sections of vibrational (de-)excitation from the vibrational level (110) to several other levels (υ′1υ′2υ′3) on the left. The corresponding vibrational levels are depicted with the same colour code on the right.

Fig. 7 Comparison of cross sections available in the literature for the excitation of the bending (010) ← (000) (top panel) and stretching (100) + (001) ← (000) (low panel) modes. The black solid line is the theoretical result obtained in this work and the blue line is from the previous study Ayouz et al. (2021).

Fig. 8 Probabilities for several rovibrational (de-)excitation of H2O. The colour code is the same as in Figs. 5 and 6. Circles of the same colour refer to the excitation process while solid or dashed lines the opposite. Note: Peff (v′ ← v) is only weakly dependent on x (i.e. temperature). It is used to obtain a cubic polynomial fit P(v' ← v) in Eq. (10) for each transition v' ← v. For convenience, the vertical dashed line at x = 0 gives the coefficient value of a0 at 300 K, evaluated by the analytical fit as ${10}^{a_0}$; see Eq. (8).

3.2 Rate coefficients

Thermally averaged rate coefficient α_υ′←υ(T) is obtained from the cross section as:

$${\alpha }_{v^{\prime }\leftarrow v}(T)=\sqrt{\frac{8}{\pi m{\left(k_{\text{b}}T\right)}^3}}{\displaystyle {\int }_0^{\infty }{\sigma }_{v^{\prime }\leftarrow v}}\left(E_{\text{el}}\right)\exp \left(-\frac{E_{\text{el}}}{k_{\text{b}}T}\right)E_{\text{el}}\text{d}{E}_{\text{el}},$$(7)

where k_b and T are the Boltzmann coefficient and the temperature, respectively.

For simplicity, we fit the obtained rate coefficients to an analytical representation in the following way. First, the quantity:

$$P^{\text{eff }}\left(v^{\prime }\leftarrow v\right)={\log }_{10}\left[{\alpha }_{v^{\prime }\leftarrow v}(T)\sqrt{T}\exp \left(\frac{{\Delta }_{v^{\prime },v}}{T}\right)\right]$$(8)

is introduced for each transition υ′←υ. It represents roughly the logarithm of overall probability for the (de-)excitation process. It is a smooth function of temperature and displayed in Fig. 8 for a few pairs of initial υ = (000) and final υ' vibrational levels as a function of x = log₁₀ (T/300 K). Thus, similarly to the method used in previous studies Kokoouline et al. (2010); Ayouz & Kokoouline (2016, 2019); Jiang et al. (2019), the numerical rate coefficients can be obtained using the following analytical formula:

$${\alpha }_{v^{\prime }\leftarrow v}^{\text{fit}}(T)=\frac{1}{\sqrt{T}}\exp \left(-\frac{{\Delta }_{v^{\prime },v}}{T}\right){10}^{P\left(v^{\prime }\leftarrow v\right)},$$(9)

where

$$P\left(v^{\prime }\leftarrow v\right)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3.$$(10)

Here, the coefficients a_i (i = 0, 1, 2, 3) are the fitting parameters. The quantity 10^{P(υ′←υ)} is roughly the (de-)excitation probability for transition υ′←υ. In Eq. (9), Δ_υ′,υ is the threshold energy defined as

$${\text{Δ}}_{v^{\prime },v}=\{\begin{array}{cc}{E}_{v^{\prime }}-E_v>0& \text{ for excitation, }\\ 0& \text{ for de-excitation. }\end{array}$$(11)

To use the formula of Eq. (9) for each pair of transitions (υ′) ↔ (υ), the temperature, T, should be in kelvins and the numerical parameters of a_i are listed in Tables A.1–A.10. For example, at T = 300 K, the (de-)excitation probability is readily given by ${10}^{a_0}$ (because x = 0 for T=300 K).

The accuracy of the fits is illustrated in Fig. 9, while Fig. 10 shows a sample of rate coefficients for the initial state (000). As observed for the cross sections, the largest rate coefficients are those for one-quantum (de)excitations.

Figure 11 compares the present rate coefficients with those by Faure and Josselin (Faure & Josselin 2008) for four low-lying transitions, which have been widely employed in astrophysical models (see e.g. Gray et al. 2016). The data from Faure & Josselin (2008), covering five vibrational levels (000), (010), (020), (100), and (001), are based on the calculations by Nishimura and Gianturco Nishimura & Gianturco (2004) for the transitions from (000) to (010), (100), and (001). The rate coefficient for the transitions (000) → (020) was derived from the theoretical rate coefficient for (000) → (010) and the propensity rule extracted from an experiment on the vibrational relaxation of H₂0 by H₂ Zittel & Masturzo (1991), as detailed in Faure & Josselin (2008). It can be observed that for one-quantum transitions the agreement between the two sets is reasonable and within a factor of 3, or better. On the other hand, for the two-quanta transition (000) → (020), the rate coefficients differ by almost two orders of magnitude. This likely reflects the different interaction types between H₂O-electron and H₂O−H₂. We conclude that specific calculations for individual state-to-state vibrational transitions are crucial and that extrapolation from different collider types is hazardous.

Fig. 9 Quality of the fit of Eq. (10). Numerical and fitted values are displayed by solid and dashed lines according the colour code of Figs. 5 and 6, respectively.

Fig. 10 Rate coefficients for several transitions starting from the ground vibrational level.

Fig. 11 Rate coefficients for the excitation of low-lying vibrational transitions. The present results (solid lines) are compared to those of Faure & Josselin, shown as dashed lines (Faure & Josselin 2008).

4 Uncertainty estimations

In this study, the main source of uncertainty is due to electron scattering calculations. To assess the uncertainty of the obtained cross sections and thermal rate coefficients, we performed scattering calculations while varying the parameters of the model. The main scattering model (model 1) is described above. In a second set of calculations (model 2), a different (smaller) basis set DZP (double-zeta polarized) and the same active space, as in model 1, were used. In model 3, the 1a, 2a, and 3a molecular orbitals were frozen, which resulted in a reduced complete active space in the configuration interaction calculations with respect to model 1 by three orbitals. A comparison of cross sections obtained using the three models is depicted in Fig. 12. Model 1 is expected to be the most accurate, because a larger basis and a larger active space were employed. From the figure, the uncertainty in the cross section for the excitation of the bending mode is changing from about 30% at energies 0.2 eV to about 10% at and above 2 eV. For the stretching modes, the uncertainty is smaller, changing from approximately 20% to 7% in the same energy interval.

Fig. 12 Uncertainties on the present theoretical results, when comparing the cross-sections for the excitation of the bending (top panel) and the stretching (bottom panel) modes, obtained using the three different models, as described in the text.

5 Conclusions

The main results of the present study are as follows:

• Cross sections for the electron-impact excitation of the 13 lowest vibrational levels of the water molecule were computed using an ab initio approach. The highest considered vibrational level was (101), located at 7200 cm⁻¹ above ground level.

• Thermally averaged rate coefficients were derived from the calculated cross sections for temperatures in the 10– 10 000 K interval and analytical fits for rate coefficients were provided.

• Uncertainty estimations of the obtained data were performed. The obtained data and their uncertainties are important for such future studies; for example, in modelling the non-LTE spectra of water in various astrophysical environments.

• Over the course of this study, we also developed a theoretical approach that can be later applied in studies of the electron-molecule excitation of other triatomic molecules and molecular ions.

Finally, we note that rotationally resolved vibrational rate coefficients are needed to run these astrophysical models Gray et al. (2016); Faure & Josselin (2008). They are currently being computed and will be reported in a future work.

Acknowledgements

This work was supported by the National Science Foundation, Grant No. PHY-2102188. It has also received funding from the program “Accueil des chercheurs étrangers” of CentraleSupélec and Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU. AF acknowledges support by the French Agence Nationale de la Recherche (ANR-Waterstars), grant number ANR-20-CE31-0011.

Appendix A Fitting parameters for the computed rate coefficients

The tables below provide fitting parameters for the computed rate coefficients.

References

All Tables

All Figures

	Fig. 1 Employed normal coordinate grid (top panel) to evaluate the integral in Eq. (3). In the lower panel, the corresponding bond coordinates grid is represented.
In the text

Fig. 2 Minimum energy path of H2O (lower panel) in eV as a function of r1 and the second bond length r2 and the bond angle θ (top panel). The bond length, r2, is given with the solid line and with values indicated on the left axis, while the bond angle is given with the dashed line with values on the right axis of the top panel. The vertical dashed line represents the position of the global minimum.

In the text

	Fig. 3 Wave functions of several lowest vibrational levels (below 7000 cm−1) of H2O. Wave functions are represented as intensity plots depending on the two (out of the three) coordinates while the third one fixed at the value 0.8. On the top of each panel, the vibrational levels and normal coordinates are indicated.
In the text

	Fig. 4 Figure details similar to Fig. 3, with several other wave functions.
In the text

	Fig. 5 Cross sections for vibrational excitation from the vibrational level (000) to several other levels (υ′1υ′2υ′3) shown on the left. H2O low-lying vibrational levels are depicted with the same colour code on the right.
In the text

	Fig. 6 Cross sections of vibrational (de-)excitation from the vibrational level (110) to several other levels (υ′1υ′2υ′3) on the left. The corresponding vibrational levels are depicted with the same colour code on the right.
In the text

	Fig. 7 Comparison of cross sections available in the literature for the excitation of the bending (010) ← (000) (top panel) and stretching (100) + (001) ← (000) (low panel) modes. The black solid line is the theoretical result obtained in this work and the blue line is from the previous study Ayouz et al. (2021).
In the text

Fig. 8 Probabilities for several rovibrational (de-)excitation of H2O. The colour code is the same as in Figs. 5 and 6. Circles of the same colour refer to the excitation process while solid or dashed lines the opposite. Note: Peff (v′ ← v) is only weakly dependent on x (i.e. temperature). It is used to obtain a cubic polynomial fit P(v' ← v) in Eq. (10) for each transition v' ← v. For convenience, the vertical dashed line at x = 0 gives the coefficient value of a0 at 300 K, evaluated by the analytical fit as ${10}^{a_0}$; see Eq. (8).

In the text

	Fig. 9 Quality of the fit of Eq. (10). Numerical and fitted values are displayed by solid and dashed lines according the colour code of Figs. 5 and 6, respectively.
In the text

	Fig. 10 Rate coefficients for several transitions starting from the ground vibrational level.
In the text

	Fig. 11 Rate coefficients for the excitation of low-lying vibrational transitions. The present results (solid lines) are compared to those of Faure & Josselin, shown as dashed lines (Faure & Josselin 2008).
In the text

	Fig. 12 Uncertainties on the present theoretical results, when comparing the cross-sections for the excitation of the bending (top panel) and the stretching (bottom panel) modes, obtained using the three different models, as described in the text.
In the text