Rotational excitation of mono- and doubly-deuterated water by hydrogen molecules

Abstract

Rate coefficients for rotational transitions in HDO and D₂O induced by H₂ collisions below 300 K are presented. Calculations have been performed at the close-coupling and coupled-states levels with the deuterated variants of the H₂O–H₂ interaction potential of Valiron et al. The HDO–H₂ and D₂O–H₂ rate coefficients are compared to the corresponding rate coefficients for HDO–He and H₂O–H₂, respectively. Significant differences are observed. In particular the new HDO rate coefficients are found to be significantly larger (by up to three orders of magnitude) than the corresponding HDO–He rate coefficients. The impact of the new HDO rate coefficients is examined with the help of non-LTE radiative transfer calculations. A number of potential HDO maser lines are finally identified, in particular the 80.6 GHz (⁠⁠) transition.

molecular data, molecular processes, stars: formation, ISM: molecules

1 INTRODUCTION

Despite the low cosmic deuterium abundance, ⁠, a spectacular deuterium enrichment of various interstellar molecules is observed in star-forming regions (Ceccarelli et al. 2007, and references therein). Among these molecules, the deuterated isotopologues of water, HDO and D₂O are of special importance because they can help to understand the origin of water in the interstellar medium and its possible link with the D/H ratios observed in comets and in the Earth’s oceans (∼ 10⁻⁴) (Hartogh et al. 2011). While interstellar HDO was detected a few years after the discovery of H₂O in the interstellar medium (Turner et al. 1975), D₂O has been identified only recently towards the protostar IRAS 16293−2422 (Butner et al. 2007). Investigation of the abundance of water and its isotopologues in space is one of the main targets of the heterodyne instrument for the infrared (HIFI) aboard Herschel Space Observatory, which was launched in 2009 May. Recent HIFI detections towards star-forming regions include HDO (Ceccarelli et al. 2010a; Comito et al. 2010; Coutens et al. 2011) and D₂O (Vastel et al. 2010) as well as the first (tentative) identification of HD¹⁸O (Bergin et al. 2010).

As the densities in star-forming regions are low (n≲ 10⁸ cm⁻³), the molecular populations are generally not at local thermodynamical equilibrium (LTE) but rather they are determined by a competition between radiative and collisional processes. The computation of molecular abundances therefore relies on radiative transfer modelling which involves the knowledge of rate coefficients for collisional (de-)excitation. As it is extremely difficult to determine these coefficients experimentally, radiative transfer models can only rely on theoretical predictions. Until recently, the only available collisional rate coefficients for deuterated water were those of Green for HDO with He (Green 1989) and those of Faure et al. for electron-impact excitation of HDO and D₂O (Faure, Gorfinkiel & Tennyson 2004). However, in regions where HDO and D₂O are abundant, the main colliding partner is H₂.

Here we report the computation of collisional rate coefficients for the rotational (de-)excitation of HDO and D₂O by H₂ molecules. Preliminary calculations have been reported by Scribano, Faure & Wiesenfeld (2010) for D₂O–H₂ and by Wiesenfeld, Scribano & Faure (2011) for HDO–H₂. These authors have presented detailed comparisons between the three water–hydrogen isotopologues H₂O–H₂, HDO–H₂ and D₂O–H₂. Significant differences were observed and they were attributed to symmetry, kinematics and intramolecular geometry effects. Moreover, in the case of HDO, rate coefficients with H₂ were found significantly larger than the HDO–He rate coefficients of Green (1989) which are currently employed in astronomical models. The aim of this paper is (i) to present the full set of collisional rate coefficients for HDO and D₂O and (ii) to investigate the impact of these new rate coefficients on non-LTE radiative transfer calculations. A brief description of the scattering calculations is presented in Section 2 along with a presentation of the rate coefficients. The general impact of the new HDO rate coefficients on radiative transfer computations is presented in Section 3. The specific case of maser transitions is discussed in Section 4. Conclusions are drawn in Section 5.

2 SCATTERING CALCULATIONS

All scattering calculations were performed using the full nine-dimensional potential energy surface (PES) for H₂O–H₂ (Faure et al. 2005; Valiron et al. 2008). This ab initio PES was obtained by combining standard coupled-cluster [CCSD(T)] calculations with explicitly correlated CCSD(T)–R12 calculations. Full details can be found in Valiron et al. (2008). This PES is independent of nuclear masses and can be employed for any water–hydrogen isotopologue. Its high accuracy has been spectacularly confirmed very recently by a number of comparisons between theory and experiment including inelastic total and differential cross-sections (Yang et al. 2010, 2011), pressure broadening cross-sections (Wiesenfeld & Faure 2010), elastic integral cross-sections (Belpassi et al. 2010) and the spectrum of the van der Waals complex (van der Avoird & Nesbitt 2011; van der Avoird et al. 2011; Wang & Carrington 2011).

The HDO–H₂ and D₂O–H₂ rigid-rotor PESs were obtained from the full flexible H₂O–H₂ potential, but in the principal inertia axes of the water isotopologues and fixing the internal coordinates at their vibrationally averaged values. Full details can be found in Scribano, Faure & Wiesenfeld (2010) and Wiesenfeld, Scribano & Faure (2011).

The quantum theory for scattering of an asymmetric top with a linear molecule can be found in Phillips, Maluendes & Green (1995). The rotational energy levels of the water isotopologues are labelled by three numbers: the angular momentum j₁ and the pseudo-quantum numbers k_a and k_c, which correspond to the projection of j₁ along the principal inertia a- and c-axes. The rotational levels of the projectile H₂ are labelled by the angular momentum j₂. In the present work, calculations were performed using the OpenMP version of the scattering code molscat.1 All calculations for D₂O were performed at the close-coupling (CC) level while for HDO the coupled-states (CS) approximation, which approximates certain Coriolis-coupling terms, was employed at the highest energies (see below). In the case of D₂O, the initial H₂ was restricted to its lowest para state (j₂= 0) because the fraction of ortho-H₂ (j₂= 1) is expected to be negligible in the coldest interstellar regions (T < 30 K) where D₂O is abundant (Troscompt et al. 2009; Vastel et al. 2010). In contrast, HDO is observed both in cold and warm (T > 100 K) environments (e.g. Coutens et al., in preparation). In these latter cases, the fractions of ortho-H₂ in j₂= 1 and para-H₂ in j₂= 2 are expected to be significant. In the case of HDO, the initial H₂ was therefore taken in its lowest para (j₂= 0, 2) and ortho (j₂= 1) levels.

All calculations included several energetically inaccessible (closed) channels to ensure that cross-sections were converged to within 5 per cent for all transitions involving levels with energies below ∼100 cm⁻¹. For higher levels in HDO, cross-sections were converged to within 20 per cent. The lowest 12 D₂O levels were considered, i.e. up to the para level 3_1,2 at 89.97 cm⁻¹ and up to the ortho level 3_1,3 at 74.51 cm⁻¹. We note that since the ortho and para levels of D₂O do not interconvert in inelastic collisions, they were treated separately. The basis set for D₂O incorporated all target states with j₁≤ 6. In the case of HDO, the lowest 30 levels were considered, i.e. up to the level 6_1,6 at 308.62 cm⁻¹. The basis set for HDO incorporated all target states up to j₁≤ 9 (limited to an energy of 1100 cm⁻¹) for CC calculations and up to j₁≤ 11 (limited to an energy of 1710 cm⁻¹) for CS calculations. The inclusion of the j₂= 2 level of H₂, which opens at 365 cm⁻¹, was found to be necessary at all collision energies for both isotopologues. In contrast, the level j₂= 3, which opens at 730 cm⁻¹, was found to have a minor influence (≤10 per cent) and was therefore neglected.

CC calculations for D₂O were performed at total energies E_tot≤ 600 cm⁻¹. In the case of HDO, CC calculations were performed at E_tot≤ 338 cm⁻¹ for para-H₂ (j₂= 0, 2) and E_tot≤ 551 cm⁻¹ for ortho-H₂ (j₂= 1). At higher energies, the CC method becomes prohibitively expensive and the CS approximation was adopted up to E_tot= 15000 cm⁻¹. The CS cross-sections were however scaled to match the CC results at the highest value of the overlapping energy range, as in Wiesenfeld et al. (2011). We also note that a few number of transitions have zero cross-sections within the CS approximation. For these, the CC results were simply extrapolated by a constant. Full details can be found in Scribano et al. (2010) and Wiesenfeld et al. (2011). Rate coefficients were obtained by integrating the cross-sections over Maxwell–Boltzmann distributions of relative velocities. For HDO colliding with para-H₂, the rate coefficients were also summed at each temperature over the final levels of H₂ (⁠⁠) and averaged over the initial levels of H₂ (j₂= 0, 2), assuming these latter have a thermal distribution. In practice, the contribution of the initial level j₂= 2 was found to be significant when the temperature exceeds about 100 K.

The rate coefficients for D₂O–H₂ and HDO–H₂ were obtained in the temperature ranges 5–100 and 5–300 K, respectively. The typical temperature dependences of the rate coefficients, the discussion of propensity rules and the comparison between para-H₂ (j₂= 0) and ortho-H₂ (j₂= 1) can be found in Wiesenfeld et al. (2011). We note in particular that in the case of HDO, the shift of the centre of mass is crucial because it breaks the C_2v symmetry and induces additional (a-type) transitions with significant rate coefficients (Wiesenfeld et al. 2011). An overall presentation of the D₂O rate coefficients at 20 K is given in Fig. 1. In this plot, the D₂O rate coefficients are plotted versus the corresponding H₂O rate coefficients for de-excitation transitions involving the lowest 12 rotational levels. The H₂O rate coefficients with para-H₂ (j₂= 0) were taken from Dubernet et al. (2006, 2009). It should be noted that these H₂O rate coefficients were obtained at the same level of accuracy as the present calculations using the H₂O–H₂ PES of Faure et al. (2005). We also note that the ortho and para symmetries are reversed in D₂O with respect to H₂O so that transitions in ortho-D₂O are compared to the corresponding ones in para-H₂O and vice versa. It can be observed that D₂O and H₂O rate coefficients generally agree within a factor of 2–3, but differences as large as a factor of ∼ 8 exist. These differences are caused by the isotopic substitution only and they were attributed by Scribano et al. (2010) to both kinematics (i.e. mass and velocities) and intramolecular geometry effects. We conclude that H/D isotopic substitution has a significant effect (for molecules of the size of H₂O) and that the rate coefficients for the main isotopologue provide only estimates for the deuterated species, at the order-of-magnitude level.

Figure 1

Comparison between the rate coefficients for D₂O and H₂O colliding with para-H₂ (j₂= 0). De-excitation transition rate coefficients are reported for the lowest 12 rotational levels at 20 K. The horizontal axis represents the present D₂O rate coefficients (in units of cm⁠) and the vertical axis represents the corresponding H₂O rate coefficients from Dubernet et al. The two solid lines delimit the region where the rate coefficients differ by less than a factor of 10. See text for details.

An overall presentation of the HDO rate coefficients at 50 and 100 K is given in Fig. 2. In this plot, the new HDO rate coefficients are plotted versus the HDO–He rate coefficients of Green (1989) for all de-excitation transitions involving the lowest 30 levels. It should be noted that the HDO–He rate coefficients were scaled by the factor of 1.35 to account for the ratio of the reduced masses between HDO–H₂ and HDO–He. A large dispersion of the rate coefficients is observed and we can note that a significant fraction of the (scaled) He rate coefficients differ from the H₂ rate coefficients by more than a factor of 10. In fact, most of the ortho-H₂ rate coefficients exceed the corresponding He rate coefficients by more than a factor of 10 and some of them by up to three orders of magnitude. We also observe that, on average, rate coefficients with ortho-H₂ (j₂= 1) are a factor of ∼5 larger than the corresponding rate coefficients for the thermalized para-H₂ (j₂= 0, 2). We conclude that He is not a good substitute for H₂ in the case of HDO, as already observed for H₂O (Faure et al. 2007; Dubernet et al. 2009). It should be noted, however, that the large differences observed here between He and H₂ partly reflect differences in the accuracy of the corresponding PES. Indeed, in the case of H₂O–He, Yang & Stancil (2007) have found that the state-to-state rate coefficients based on modern accurate PES are larger than the corresponding rate coefficients obtained from the PES of Palma et al. (1988), i.e. the one employed by Green (1989) for HDO–He, by up to a factor of 2 at 50 K. In any case, a significant impact of the new HDO rate coefficients is thus expected in the determination of interstellar HDO abundances, as examined below.

Figure 2

Comparison between the rate coefficients for HDO–H₂ and HDO–He. De-excitation transition rate coefficients are reported for the lowest 30 rotational levels at 50 and 100 K. The horizontal axis represents the present HDO–H₂ rate coefficients (in units of cm³ s⁻¹) and the vertical axis represents the corresponding (scaled) HDO–He rate coefficients from Green (1989). Left- and right-hand panels correspond to thermalized para-H₂ (j₂= 0, 2) and ortho-H₂ (j₂= 1), respectively. The two solid lines delimit the region where the rate coefficients differ by less than a factor of 10. See text for details.

We note finally that the full sets of HDO and D₂O rate coefficients are made available in the LAMDA data base2 as well as at the Centre de données astronomiques de Strasbourg (CDS).3 As explained above, the rate coefficients for HDO colliding with para-H₂ are summed and averaged assuming a thermal distribution of j₂= 0, 2. However, for applications where the H₂ rotational populations are not thermalized, the j₂-resolved rate coefficients can be obtained upon request to the authors.

3 IMPACT OF THE NEW HDO COLLISION RATE COEFFICIENTS

In order to estimate the impact of the new rate coefficients on the HDO abundance derived from observations, we have performed non-LTE radiative transfer calculations with the two sets of rate coefficients, the old ones from Green (1989) and the new ones from this work. The large velocity gradient (LVG) approximation was employed with a code adapted from the one described in Ceccarelli et al. (2003). It refers to a semi-infinite isodense and isothermal slab in plane-parallel geometry. The first 30 levels of HDO were included. The ortho-to-para ratio of H₂ is taken at its thermodynamic equilibrium value for each temperature. As a background radiation field, we included the 2.7 K cosmic microwave background only. Fig. 3 gives an overview of the results for all lines in the 0–2000 GHz frequency range and with an upper level energy lower than 150 K. Specifically, we plot the ratio between the brightness temperature T_B values calculated with the new rate coefficients and those calculated with the Green (1989) rate coefficients, for different densities and temperatures, in the optically thin case. Fig. 3 makes it very clear that this ratio is between 3 and 30 for most of the quoted transitions and that it can be as high as 100 for a few of them, in the range of temperature (30–200 K), densities (10⁵–10⁸ cm⁻³) and HDO column density (10¹⁴–10¹⁷ cm⁻²) explored here. We conclude that in the optically thin case, the inaccuracies in excitation rate coefficients are quasi-linearly propagated within the radiative transfer equations. This implies a possible overestimation of the HDO abundance (and therefore of the D/H ratio of water) in the emitting region by the same value, depending on the used line.

Figure 3

Ratio of the brightness temperature T_B values of different transitions calculated with the present rate coefficients with respect to the Green (1989) rate coefficients, as a function of the upper level energy of the transition. Different colours refer to different temperatures: 30 K (red), 100 K (green) and 200 K (blue). Different symbols refer to different densities: 10⁵ cm⁻³ (diamonds), 10⁶ cm⁻³ (triangles), 10⁷ cm⁻³ (squares) and 10⁸ cm⁻³ (crosses). In these computations, a HDO column density of 10¹⁴ cm⁻² (km s⁻¹)⁻¹ was adopted.

Fig. 3 also shows that the largest ratios are obtained at the lowest H₂ densities and temperatures. This is expected since the farther from the LTE, the larger the importance of the employed excitation rate coefficients on the level population. However, this is not strictly true in all situations. Some specific cases of used HDO lines are reported in Figs 4–6, in order to show how exactly the ratio depends on the density and temperature of the gas as well as on the HDO column density (namely the line’s optical depth). Fig. 4 shows the specific case of the HDO line at 893 GHz. The impact of the new rate coefficients on this line is large at low temperatures and low densities, as the ratio can be a factor of almost 30. The largest difference occurs at the lowest HDO column density (left-hand panel), where the line optical depths do not influence the levels population. At large enough HDO column densities, the ratio tends to be less than 2. Figs 5 and 6 report two cases where the impact of the new collisional coefficients is particularly large, and the difference between the present and old collisional coefficients leads to factors of larger than 100. In the case of the 241-GHz line, often observed, the largest difference between the two sets of rate coefficients is not necessarily obtained at the lowest column densities (left-hand panel). Indeed, the maximum difference is found to peak at ∼ 2 × 10¹⁵ cm⁻². It should be noted that the oscillation of the plotted ratio in the range 10¹⁴–10¹⁵ cm⁻² is not unphysical: the brightness temperatures were checked and they show a perfectly smooth dependence with respect to the column density. This oscillation simply reflects the great sensitivity of radiative transfer equations to the excitation rate coefficients in this range of column density. We finally note that for the majority of the lines, the largest difference between the two sets of rate coefficients peaks at densities around 10⁵ cm⁻³, which is consistent with the typical value of the critical density, n_cr∼ 10⁷ cm⁻³.

Figure 4

Transition at frequency 893.6 GHz. Ratio of the T_B values calculated with the present rate coefficients with respect to the Green (1989) rate coefficients as a function of the HDO column density (left-hand panel) and H₂ density (right-hand panel). Different colours refer to different temperatures: 30 K (red), 100 K (green) and 200 K (blue). In the left-hand panel, different lines refer to different H₂ densities: 10⁵ cm⁻³ (solid line), 10⁶ cm⁻³ (dotted line), 10⁷ cm⁻³ (dashed line) and 10⁸ cm⁻³ (dash–dotted line). In the right-hand panel, different lines refer to different HDO column densities: 10⁴ cm⁻² (solid line), 10¹⁵ cm⁻² (dotted line), 10¹⁶ cm⁻² (dashed line) and 10¹⁷ cm⁻² (dash–dotted line).

Figure 5

Same as Fig. 4, but for the transition at frequency 599.9 GHz.

Figure 6

Same as Fig. 4, but for the transition at frequency 241.6 GHz.

4 MASER EMISSION

In this section, we discuss the specific cases of the lines which present a potential maser effect, using the same LVG code of Section 3. Transitions with inverted level population (masing lines) are characterized by negative optical depths and negative excitation temperatures and the strongest masers have the largest negative optical depths. Fig. 7 shows the optical depth τ of the transition at 225.9 GHz. Negative values of τ are observed at 100 and 200 K, for column densities above 10¹⁴ cm⁻². The largest negative optical depth, τ=− 1.1, is reached for column densities between 2 and 8 × 10¹⁶ cm⁻² and a H₂ density between 2 and 20 × 10⁶ cm⁻³. The corresponding brightness temperatures, T_B, are plotted in Fig. 8. It is observed that T_B can reach 1200 K for a HDO column density of 10¹⁷ cm⁻² and H₂ densities between 2 and 8 × 10⁶ cm⁻³.

Figure 7

Optical depth τ of the transition at 225.9 GHz. Different colours refer to different temperatures: 30 K (red), 100 K (green) and 200 K (blue). In the left-hand panel, different lines refer to different H₂ densities: 10⁵ cm⁻³ (solid lines), 10⁶ cm⁻³ (dotted lines), 10⁷ cm⁻³ (dashed lines) and 10⁸ cm⁻³ (dash–dotted lines). In the right-hand panel, different lines refer to different HDO column densities: 10¹⁴ cm⁻² (solid lines), 10¹⁵ cm⁻² (dotted lines), 10¹⁶ cm⁻² (dashed lines) and 10¹⁷ cm⁻² (dash–dotted lines).

Figure 8

Brightness temperature T_B (in K) for the transition at 225.9 GHz. Lines and colours are as in Fig. 7.

Six other lines were found to invert in the frequency range 0–2000 GHz, notably the often observed 80.6- and 255.1-GHz lines. They are reported in Table 1. The inversions are found to occur for similar physical conditions: T∼ 100–200 K,  cm⁻² and  cm⁻³. The maximum brightness temperatures, T_B= 1200 K, are observed for the two transitions at 80.6 and 225.9 GHz. We note that while the 225.9-GHz line also inverts with the excitation rate coefficients of Green (1989) (but at larger densities), the 80.6-GHz line, as well as the 5.7-GHz line, was not predicted to invert with the collisional coefficients of Green (1989). This is a particularly interesting result, as the 80.6-GHz line is observable and observed from ground telescopes. The masering effect, not previously identified, can help in explaining the observations obtained by Codella et al. (2010) towards the low-mass protostar L1448-mm. They found the 80.6-GHz emission associated with the hot corino of the source and a spot possibly associated with the outflow shock. Intriguingly, they noticed that the spatial position of this spot coincided with the distribution of the H₂O 22 GHz and, possibly, 183-GHz maser lines. The new collisional coefficients support the interpretation that the HDO observed emission is associated with regions of high densities, where the line mases.

Table 1

HDO transitions presenting a population inversion. The ranges of column density, temperature and H₂ density where the optical depth τ is negative is given for each transition. The largest values of |τ | and T_B (brightness temperature) are also provided.

Frequency (GHz)	Transition	Density (cm−3)	Temperature (K)	N (HDO) (cm−2)	|τmax |	(K)
0.8			100–200	> 2 × 1016	3	60
5.7			100–200	> 1016	1	18
80.6			100–200	> 2 × 1015	1	1200
120.8			100–200	> 2 × 1015	2.5	550
138.5			100–200	> 2 × 1016	2	150
225.9			100–200	> 1015	1	1200
255.1			100–200	> 1016	1.5	450

As a final remark, we note that the conditions for the HDO lines to mase (relatively high densities, temperature and HDO column densities) can be found in several astronomical situations. First, in the hot cores and hot corinos of protostars, but also in the outflow shocks and the protoplanetary discs’ innermost regions. Indeed, Ceccarelli et al. (2010b) noted that the 255.1-GHz line mases, by using the Green (1989) collisional coefficients, in protoplanetary discs.

5 CONCLUSIONS

We have reported rate coefficients for rotational transitions in HDO and D₂O induced by H₂ collisions below 300 K. Calculations were performed at the CC and CS levels with the deuterated variants of the H₂O–H₂ interaction potential of Valiron et al. (2008). The HDO–H₂ and D₂O–H₂ rate coefficients were compared to the corresponding rate coefficients for HDO–He (Green 1989) and H₂O–H₂ (Dubernet et al. 2006, 2009), respectively. Significant differences were observed and it was concluded (i) that the rate coefficients for H₂O provide only estimates for D₂O at the order-of-magnitude level and (ii) that He is not a good substitute for H₂ in the case of HDO, with differences up to three orders of magnitude. The impact of the new HDO rate coefficients was examined using non-LTE LVG modelling and it was shown that in the optically thin case the inaccuracies in excitation rate coefficients are quasi-linearly propagated within the radiative transfer equations. The present rate coefficients should therefore lead to a significant re-estimation of the HDO abundance and, therefore, of the D/H ratio of interstellar water (see e.g. Coutens et al., in preparation). A number of potential maser transitions were also identified in the frequency range 0–300 GHz, with optical depths larger than 1 in absolute value and brightness temperatures up to 1200 K. A particularly interesting result is the prediction of a masering effect in the 80.6-GHz line, supporting the interpretation that the HDO observed emission towards the low-mass protostar L1448-mm is associated with regions of high densities (Codella et al. 2010).

Finally, we note that in contrast to H/D isotopic substitution, isotopic effects for the oxygen isotopologues of water (HO and HO) can be safely ignored since the change in internal geometries is essentially negligible and the C_2v symmetry is conserved.

Footnotes

1

Repository at http://ipag.osug.fr/~afaure/molscat/index.html.

2

http://www.strw.leidenuniv.nl/~moldata/ (Schöier et al. 2005)

3

http://cdsweb.u-strasbg.fr

This work has been supported by the CNRS national programme ‘Physique et Chimie du Milieu Interstellaire’ (PCMI). Scattering calculations were performed on the CIMENT supercomputing platform in Grenoble with valuable help from F. Roch.