# SOFIA/GREAT observations of OD and OH rotational lines towards high-mass star forming regions

T. Csengeri, F. Wyrowski, K.M. Menten, H. Wiesemeyer, R. Güsten, J.

Stutzki, S. Heyminck, Y. Okada

► To cite this version:

T. Csengeri, F. Wyrowski, K.M. Menten, H. Wiesemeyer, R. Güsten, et al. SOFIA/GREAT observations of OD and OH rotational lines towards high-mass star forming regions. Astronomy and Astrophysics - A&A, 2022, 658, pp.A193. 10.1051/0004-6361/202140577. hal-03531114

HAL Id: hal-03531114

https://hal.science/hal-03531114

Submitted on 11 Apr 2023

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
# SOFIA/GREAT observations of OD and OH rotational lines towards high-mass star forming regions

T. Csengeri<fnref n="1" />,<fnref n="2" />, F. Wyrowski<fnref n="2" />, K. M. Menten<fnref n="2" />, H. Wiesemeyer<fnref n="2" />, R. Güsten<fnref n="2" />, J. Stutzki<fnref n="3" />, S. Heyminck<fnref n="2" />, and Y. Okada<fnref n="3" />

Received 16 February 2021 / Accepted 26 November 2021

## ABSTRACT

**Context.** Only recently, OD, the deuterated isotopolog of hydroxyl, OH, has become accessible in the interstellar medium; spectral lines from both species have been observed in the supra-Terahertz and far infrared regime. Studying variations of the OD/OH abundance amongst different types of sources can deliver key information on the formation of water, H₂O.

**Aims.** With observations of rotational lines of OD and OH towards 13 Galactic high-mass star forming regions, we aim to constrain the OD abundance and infer the deuterium fractionation of OH in their molecular envelopes. For the best studied source in our sample, G34.26+0.15, we were able to perform detailed radiative transfer modelling to investigate the OD abundance profile in its inner envelope.

**Methods.** We used the Stratospheric Observatory for Infrared Astronomy (SOFIA) to observe the <fnref n="2" />Π₃/₂ J = 5/2–3/2 ground-state transition of OD at 1.3 THz (215 µm) and the rotationally excited OH line at 1.84 THz (163 µm). We also used published high-spectral-resolution SOFIA data of the OH ground-state transition at 2.51 THz (119.3 µm).

**Results.** Absorption from the <fnref n="2" />Π₃/₂ OD J = 5/2–3/2 ground-state transition is prevalent in the dense clumps surrounding active sites of high-mass star formation. Our modelling suggests that part of the absorption arises from the denser inner parts, while the bulk of it as seen with SOFIA originates in the outer, cold layers of the envelope for which our constraints on the molecular abundance suggest a strong enhancement in deuterium fractionation. We find a weak negative correlation between the OD abundance and the bolometric luminosity to mass ratio, an evolutionary indicator, suggesting a slow decrease of OD abundance with time. A comparison with HDO shows a similarly high deuterium fractionation for the two species in the cold envelopes, which is of the order of 0.48% for the best studied source, G34.26+0.15.

**Conclusions.** Our results are consistent with chemical models that favour rapid exchange reactions to form OD in the dense cold gas. Constraints on the OD/OH ratio in the inner envelope could further elucidate the water and oxygen chemistry near young high-mass stars.

**Key words.** stars: formation – HII regions – ISM: molecules – submillimeter: ISM

## 1. Introduction

The hydroxyl radical, OH, is a key molecule to our understanding of interstellar chemistry, in particular to constrain chemical pathways leading to the formation of water, a prime molecule of interest for star- and planet formation (for a review see van Dishoeck et al. 2013). Tracing its deuterated form, OD, promises to place additional constraints on formation routes through the establishment of the deuterium fractionation within protostellar envelopes.

The OH radical has been observed in various environments, first through its radio lines corresponding to transitions between rotational energy levels split by the magnetic hyperfine structure (hfs) of the molecule (see e.g. Weinreb et al. 1963 and Weaver & Williams 1964). Cesaroni & Walmsley (1991) discussed observations of multiple hfs transitions and modelled their excitation, which generally is not governed by local thermodynamic equilibrium (LTE). As shown in their energy level diagrams in Fig. 1, the rotational transitions of OH (and those of OD) lie in the supra-THz regime at >1.3 THz and are only accessible from space and with airborne facilities. The first observations

of rotational transitions of OH were performed with the *Kuiper Airborne Observatory* (Storey et al. 1981; Melnick et al. 1987; Betz & Boreiko 1989) and have been followed by the Infrared Space Observatory (Ceccarelli et al. 1998; Giannini et al. 2001; Nisini et al. 2002), all with modest spectral resolution of tens of km s⁻<fnref n="1" />. Only recently, observations of rotational transitions of OH have become accessible at higher spectral resolution with the Heterodyne Instrument for the Far-Infrared aboard *Herschel* (HIFI; Wampfler et al. 2011, 2013) and the German Receiver for Astronomy at Terahertz Frequencies (GREAT) on SOFIA (Csengeri et al. 2012; Wiesemeyer et al. 2012, 2016; Parise et al. 2012). With this receiver, the ground-state rotational line of its deuterated form, OD, was first detected in absorption towards the envelope of the low-mass protostar IRAS 16293–2422 (Parise et al. 2012).

The OH molecule is a key component in the water chemistry pathway. At low temperatures, both species are intimately related, forming in the gas phase via dissociative recombination of the oxonium ion, H₃O⁺ (Neau et al. 2000; Jensen et al. 2000 and references therein). The latter is synthesised through a sequence of molecular hydrogen abstraction reactions

<footnote n="1">¹ Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France
e-mail: timea.csengeri@u-bordeaux.fr</footnote>

<footnote n="2">² Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote n="3">³ I. Physikalisches Institut der Universität zu Köln, Zülpicher Strasse 77, 50937 Köln, Germany</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Ground state and first excited energy levels of the $^2\Pi_{3/2}$ OD, OH, and $^2\Pi_{1/2}$ OH states. Grey arrows mark the transitions discussed in this study, whose frequencies are given at the bottom. The zero level of the energy scale refers to the vibrational ground states (2658 and 1941 K for OH and OD, respectively).</figcaption></figure>

(Dalgarno 1994) and is therefore sensitive to the molecular gas fraction. At high temperatures, OH forms through gas phase reactions from the reaction between atomic oxygen and hydrogen molecules, and through successive hydrogenation is converted into water (Charnley 1997; Harada et al. 2010; van Dishoeck et al. 2014). The reaction of OH + H₂ ⇌ H₂O + H, which is endothermic by 920 K, is sensitive to the ratio of atomic O and H and thus depends strongly on the local conditions, such as the temperature and UV field. OH is also a by-product of H₂O photodissociation, making it a sensitive tracer of the radiation field. This, for example, results in the high OH abundances in the dense molecular envelopes of (ultra)compact HII regions (see e.g. Cesaroni & Walmsley 1991; Csengeri et al. 2012). The formation of OH is also dependent on the cosmic ray ionisation rate and the gas density (Wakelam et al. 2010).

Since OH is chemically linked to the formation of water, to measure the OD/OH ratio in star forming regions presents a promising tool to constrain the formation conditions of water. For the first time, we can observationally constrain this ratio, which is expected to probe the chemical formation pathway for these molecules, as deuterium enhancement (i.e. fractionation) is strongly dependent on the chemical route and its formation conditions (e.g. Croswell & Dalgarno 1985; Roberts & Herbst 2002). Since deuterium is of solely primordial origin (Peebles 1966) and is destroyed in stars through fusion to $^3$He, a better understanding of its fractionation into OD may ultimately help to set constraints on star formation history.

Making use of the capabilities of the SOFIA/GREAT receiver, we investigated the OD/OH ratio towards a sample of 13 Galactic high-mass, star forming regions. Most of our targets are well-known Galactic massive star-forming regions subject to extensive studies (e.g. van der Tak et al. 2006; Wyrowski et al. 2016; König et al. 2017). They cover a range of evolutionary stages, characterized by clumps hosting embedded ultra-compact HII (UCHII) regions and internally heated hot molecular cores, as well as quiescent massive sources. Our classification of the evolutionary stage of the sample is based on the type of embedded sources associated with them. In short, sources associated with compact radio emission from the CORNISH catalogue Purcell et al. (2013) or Walsh et al. (1998) are listed as UCHII

regions. Parallel or already prior to this stage clumps appear as infrared-bright luminous sources due to the internally heated gas, often referred to as infrared bright massive clumps (cf. Csengeri et al. 2016; König et al. 2017). This stage is frequently characterised by rich molecular emission mainly originating from complex organic molecules leading to the emergence of hot cores. Many of our targets are well-characterised rich molecular hot cores that we also indicated in Table 1. The two quiescent sources G23.21–0.3 and G34.41+0.2 lack bright emission at mid-infrared wavelengths and are likely to correspond to the earliest stages in the evolutionary sequence leading to high-mass star and cluster formation. For this sample of sources, this classification approach is in agreement with a classification based on the physical properties of the clumps, using their bolometric luminosity ($L_{bol}$) and clump mass (cf. König et al. 2017). We list the sources and give a summary of the observed lines and transitions in Table 1.

This paper is organised as follows. In Sect. 2, we introduce the observed data and the used datasets, and in Sect. 3 we present our observational results. We perform a detailed non-local thermodynamic equilibrium (non-LTE) radiative transfer modelling in Sect. 4 and discuss our results in Sect. 5.

## 2. Observations and data reduction

### 2.1. SOFIA/GREAT observations

We used the high spectral resolution GREAT receiver<fnref n="1" /> (Heyminck et al. 2012; Risacher et al. 2016) on board the Stratospheric Observatory for Infrared Astronomy (SOFIA; Young et al. 2012) to observe rotational transitions of the OD and OH lines. Our primary focus was the $^2\Pi_{3/2}$ ground-state ($J = 5/2-3/2$) transition of 2OD at 1391.5 GHz, which was complemented towards the majority of our sample by observations of the $^2\Pi_{1/2}$ excited-state ($J = 3/2-1/2$) Λ-doublet transitions of OH at 1834.7 (163.4 µm) and 1837.8 GHz (163.1 µm), respectively.

The observations were carried out as part of GREAT consortium guaranteed observing time in the observatory’s Cycle 1, and Cycle 4 programmes. Best spectroscopic baselines were obtained by means of the chop-and-nod procedure (double-beam switching), with chop amplitudes of 30′′–60′′, operated at a frequency of 1 Hz. Owing to the double-sideband reception of GREAT, the selection of the sideband and intermediate frequency was optimized with regard to the atmospheric transmission (so as to avoid ozone features) and to receiver sensitivity. In Cycle 1 for OD, the signal was recorded at 1391.5 GHz using the L1 band of GREAT tuned in the upper side-band, and for the $^2\Pi_{1/2}$ excited-state OH line using the L2 band at 1837.8 GHz tuned in the upper side-band. The 2.5 GHz wide intermediate frequency band was tuned to 1.7 GHz offset, allowing us to simultaneously cover the 1834.7 GHz Λ-doublet pair in the image band. In Cycle 4, the $^2\Pi_{1/2}$ excited-state OH line was observed using the GREAT low frequency array module (LFA) array and was also tuned to 1837.8 GHz in the upper side-band. We also recorded the $^2\Pi_{3/2}$ ground-state $^{18}$OH line at 2494.7 GHz towards G34.26+0.15 and G351.58–0.4 using the $M_\alpha$ channel of GREAT (Heyminck et al. 2012) in Cycle 1 in July 2013. The receiver was tuned to 2494.420 GHz in the upper side band.

<footnote n="1">¹ GREAT is a development by the MPI für Radioastronomie and KOSMA/Universität zu Köln, in cooperation with the MPI für Sonnensystemforschung and the DLR Institut für Optische Sensorsysteme.</footnote>
<table><caption>**Table 1.** Target sources and observed species.</caption><thead><tr><th>Source</th><th>α<br>(J2000)</th><th>δ<br>(J2000)</th><th>ν<sub>lsr</sub><br>(km s<sup>-1</sup>)</th><th>Obs. species</th><th>OD<br>detection</th><th>Type</th></tr></thead><tbody><tr><td>G10.47+0.03</td><td>18:08:38.20</td><td>-19:51:50.0</td><td>+66.8 <sup>(a)</sup></td><td>OD, OH</td><td>abs</td><td>UCHII, hot core <sup>(1,2,3)</sup></td></tr><tr><td>G23.21-0.3</td><td>18:34:54.91</td><td>-08:49:19.2</td><td>+76.8 <sup>(a)</sup></td><td>OD</td><td>abs</td><td>quiescent</td></tr><tr><td>G34.26+0.15</td><td>18:53:18.49</td><td>01:14:58.7</td><td>+58.1 <sup>(a)</sup></td><td>OD, OH, <sup>18</sup>OH</td><td>abs</td><td>UCHII, hot core <sup>(1,3)</sup></td></tr><tr><td>G34.41+0.2</td><td>18:53:18.13</td><td>01:25:23.7</td><td>+57.6 <sup>(a)</sup></td><td>OD</td><td>–</td><td>Quiescent</td></tr><tr><td>G35.20-0.7</td><td>18:58:12.93</td><td>01:40:40.6</td><td>+33.4 <sup>(a)</sup></td><td>OD, OH<sup>†</sup></td><td>abs</td><td>Hot core <sup>(4)</sup></td></tr><tr><td>W28A (G5.89-0.4)</td><td>18:00:30.40</td><td>-24:04:00.0</td><td>+11.0 <sup>(a)</sup></td><td>OD, OH</td><td>Uncertain</td><td>UCHII, hot core <sup>(3)</sup></td></tr><tr><td>W31C (G10.62)</td><td>18:10:28.69</td><td>-19:55:50.0</td><td>-2.0 <sup>(b)</sup></td><td>OD, OH</td><td>abs</td><td>UCHII, hot core <sup>(1)</sup></td></tr><tr><td>W33A</td><td>18:14:39.52</td><td>-17:52:00.5</td><td>+36.7 <sup>(a)</sup></td><td>OD, OH <sup>(†)</sup></td><td>–</td><td>UCHII, hot core <sup>(5,6)</sup></td></tr><tr><td>W49N</td><td>19:10:13.20</td><td>09:06:12.0</td><td>+12.8 <sup>(a)</sup></td><td>OD, OH</td><td>Two comp. abs</td><td>UCHII, hot core <sup>(1,3)</sup></td></tr><tr><td>W51e2</td><td>19:23:43.96</td><td>14:30:34.6</td><td>+58.2 <sup>(a)</sup></td><td>OD</td><td>Two comp. abs</td><td>UCHII, hot core <sup>(7)</sup></td></tr><tr><td>W51d</td><td>19:23:39.95</td><td>14:31:08.0</td><td>+58.2 <sup>(a)</sup></td><td>OD</td><td>–</td><td>UCHII <sup>(1,7)</sup></td></tr><tr><td>G327.29-0.6</td><td>15:53:08.55</td><td>-54:37:05.1</td><td>-44.7 <sup>(a)</sup></td><td>OD, OH</td><td>abs</td><td>Hot core <sup>(7)</sup></td></tr><tr><td>G351.58-0.4</td><td>17:25:25.03</td><td>-36:12:45.4</td><td>-95.9 <sup>(a)</sup></td><td>OD, OH, <sup>18</sup>OH</td><td>abs</td><td>UCHII, hot core <sup>(7,8)</sup></td></tr></tbody></table>

**Notes.** Columns are, left to right, source designation, J2000 coordinates, LSR velocity, observed species, information on OD detection, and information on source type. <sup>(a)</sup>Based on C<sup>17</sup>O (J = 3–2) measurements with the APEX telescope from Giannetti et al. (2014) and Wyrowski et al. (2016). <sup>(b)</sup>Based on NH<sub>3</sub> (1,1) measurements with the Effelsberg telescope from Wienen et al. (2012). <sup>(†)</sup>Only the excited OH lines have been observed.

**References.** <sup>(1)</sup> Purcell et al. (2013), <sup>(2)</sup>Rolffs et al. (2011b), <sup>(3)</sup>Hatchell et al. (1998), <sup>(4)</sup>Sánchez-Monge et al. (2013), <sup>(5)</sup>van der Tak & Menten (2005), <sup>(6)</sup>Immer et al. (2014), <sup>(7)</sup>Rolffs et al. (2011a), <sup>(8)</sup>Walsh et al. (1998).

We used the extended fast Fourier transform spectrometers (XFFTS) as backends with 2.5 GHz bandwidth and 76 kHz spectral resolution in Cycle 1, and the FFTS4G in Cycle 4 providing a 4 GHz bandwidth. The median single-sideband system temperature was ~1200 K at 1391.5 GHz, and ~1600 K at 1834.7 GHz (Rayleigh-Jeans equivalent antenna temperatures at systemic velocity). The pointing was established with the optical guide cameras providing an accuracy of 2–3′′, and linked to the instrument reference frame by peaking up on Mars or Jupiter right after each reinstallation of the GREAT receiver system. The observed lines for each source are listed in Table 1, and a detailed summary of the observed transitions is listed in Table 2.

The data were calibrated using standard procedures based on the KOSMA/GREAT calibrator (Guan et al. 2012), and for the data analysis we used the GILDAS<sup>2</sup> software. The transmission correction was based on the AM atmosphere model<sup>3</sup>, scaling the measured total power by a forward efficiency of 97%. We fine-tuned the correction for the ozone features by scaling their column densities in AM by factors of order unity until they disappeared against the continuum signal. The data was then converted to a Rayleigh-Jeans equivalent main-beam brightness temperature ($T_{mb,RJ}$) scale using a beam efficiency of 67% at both frequencies for Cycle 1, and of 70% for the 1837.8 GHz observations in Cycle 4. For the <sup>18</sup>OH line, we used a beam efficiency of 70%.

For the OD line, we obtained the continuum emission by means of a dedicated double-sideband calibration. The Λ-doublet pair of the <sup>2</sup>Π<sub>1/2</sub> excited-state OH lines originates from different sidebands and was therefore calibrated considering the different atmospheric transmission in the image and signal bands. The spectra from the signal and image band calibration were then stitched together at an emission-free velocity to obtain a composite spectrum with the respective atmospheric calibration for the spectral features from both the signal and image band.

For both the <sup>2</sup>Π<sub>3/2</sub> OD (J = 5/2–3/2) and the <sup>2</sup>Π<sub>3/2</sub> OH (J = 3/2–1/2) lines, the spectra were first summed up to define windows excluding the lines for the baseline fitting. Then we used a polynomial baseline to fit the continuum emission and subtracted the fitted continuum from each scan individually. The spectra were then averaged with noise weighting and were smoothed to ~0.65 and 0.75 km s<sup>-1</sup> at 1391.5 and 1837.8 GHz, respectively. We reach an average noise of 84 and 87 mK for these velocity resolutions. The half-power beam width (HPBW) is ~20′′ at 1391.5 GHz, 16′′ at 1837.8 GHz. The same procedure has been applied for the <sup>18</sup>OH line; however, the baseline variations needed a more careful treatment for G34.26+0.15 due to instrumental effects (see Appendix C).

## 2.2. Archival SOFIA data

To put better constraints on the origin of the emission, we used data resulting from SOFIA observations of the <sup>2</sup>Π<sub>3/2</sub>, J = 5/2–3/2 OH rotational ground-state transition at 2514 GHz from Wiesenmeyer et al. (2012, 2016). The observational setup and data reduction steps are similar to those described above. The final spectra were smoothed to 0.9 km s<sup>-1</sup> velocity resolution and have a median noise level of ~0.3 K. The HPBW is 12′′ at 2514 GHz.

## 2.3. Archival Herschel data

To compare the OD and HDO line profiles, we used spectra of the HDO (J<sub>K,A</sub> = 1<sub>1,1</sub>-0<sub>0,0</sub>) line at 893.6 GHz (Coutens et al. 2014) observed by the Herschel/HIFI instrument, where available. In addition to this, to assist the detailed modelling of G34.26+0.15, we also used a spectrum of the p-H<sup>18</sup>O (1<sub>1,1</sub>-0<sub>0,0</sub>) line at 1101.698 GHz (Flagey et al. 2013). These observations were made in the framework of the Probing Interstellar Molecules with Absorption Line Studies (PRISMAS) Guaranteed Time Key programme (Gerin et al. 2010; Neufeld et al. 2010). PRISMAS observations employed a dual-beam-switching pointed observing mode with a reference position of ~3′ from the source. We obtained the level 2.5 data products from the

<footnote n="2"><sup>2</sup> See http://www.iram.fr/IRAMFR/GILDAS/</footnote>

<footnote n="3"><sup>3</sup> S. Paine, SMA Memo #152, Smithsonian Astrophysical Observatory.</footnote>
<table><caption>**Table 2.** Molecules and their transitions.</caption><thead><tr><th>Molecule</th><th>Transition</th><th>Hyperfine component</th><th>Frequency (GHz)</th><th>Δν<sub>HFS</sub><sup>(*)</sup> (km s<sup>-1</sup>)</th><th>A<sub>E</sub><sup>(†)</sup> (s<sup>-1</sup>)</th><th>E<sub>up</sub>/k (K)</th><th>Instrument</th><th>B<sub>eff</sub></th><th>HPBW (′′)</th><th>Reference</th></tr></thead><tbody><tr><td rowspan="6">OD</td><td rowspan="6"><sup>2</sup>Π<sub>3/2</sub> J = 5/2 ← 3/2</td><td>F = 3/2<sup>-</sup> ← 5/2<sup>+</sup></td><td>1391.4857</td><td>+1.9</td><td>2.3<sup>-</sup>(-4)</td><td rowspan="6">67</td><td rowspan="6">SOFIA/GREAT</td><td rowspan="6">0.67</td><td rowspan="6">20′′</td><td rowspan="6">(1)</td></tr><tr><td>F = 5/2<sup>-</sup> ← 5/2<sup>+</sup></td><td>1391.4895</td><td>+1.1</td><td>3.7<sup>-</sup>(-3)</td></tr><tr><td>F = 7/2<sup>-</sup> ← 5/2<sup>+</sup></td><td>1391.4947</td><td>0.0</td><td>2.3<sup>-</sup>(-2)</td></tr><tr><td>F = 3/2<sup>-</sup> ← 3/2<sup>+</sup></td><td>1391.4980</td><td>-0.7</td><td>5.6<sup>-</sup>(-3)</td></tr><tr><td>F = 5/2<sup>-</sup> ← 3/2<sup>+</sup></td><td>1391.5017</td><td>-1.5</td><td>2.0<sup>-</sup>(-2)</td></tr><tr><td>F = 3/2<sup>-</sup> ← 1/2<sup>+</sup></td><td>1391.5051</td><td>-2.3</td><td>1.8<sup>-</sup>(-2)</td></tr><tr><td rowspan="4">OH</td><td rowspan="4"><sup>2</sup>Π<sub>1/2</sub> J = 3/2 → 1/2</td><td>F = 1<sup>-</sup> → 1<sup>+</sup></td><td>1834.7350</td><td>+1.9</td><td>2.1<sup>-</sup>(-2)</td><td rowspan="4">270</td><td rowspan="4">SOFIA/GREAT</td><td rowspan="4">0.67<sup>(‡)</sup></td><td rowspan="4">16′′</td><td rowspan="4">(1)</td></tr><tr><td>F = 2<sup>-</sup> → 1<sup>+</sup></td><td>1834.7469</td><td>0.0</td><td>6.4<sup>-</sup>(-2)</td></tr><tr><td>F = 1<sup>-</sup> → 0<sup>+</sup></td><td>1834.7499</td><td>-0.5</td><td>4.2<sup>-</sup>(-2)</td></tr><tr><td>F = 1<sup>+</sup> → 1<sup>-</sup></td><td>1837.7461</td><td>+11.5</td><td>2.1<sup>-</sup>(-2)</td></tr><tr><td rowspan="6"><sup>18</sup>OH</td><td rowspan="6"><sup>2</sup>Π<sub>1/2</sub> J = 3/2 → 1/2</td><td>F = 2<sup>+</sup> → 1<sup>-</sup></td><td>1837.8163</td><td>0.0</td><td>6.4<sup>-</sup>(-2)</td><td rowspan="6">270</td><td rowspan="6">SOFIA/GREAT</td><td rowspan="6">0.67<sup>(‡)</sup></td><td rowspan="6">16′′</td><td rowspan="6">(1)</td></tr><tr><td>F = 2<sup>+</sup> → 0<sup>-</sup></td><td>1837.8365</td><td>-3.3</td><td>4.3<sup>-</sup>(-2)</td></tr><tr><td>F = 1<sup>+</sup> → 0<sup>+</sup></td><td>2494.6809</td><td>+1.7</td><td>1.4<sup>-</sup>(-2)</td></tr><tr><td>F = 3<sup>-</sup> → 2<sup>+</sup></td><td>2494.6950</td><td>0.0</td><td>1.4<sup>-</sup>(-1)</td></tr><tr><td>F = 2<sup>-</sup> → 1<sup>+</sup></td><td>2494.7342</td><td>-4.7</td><td>1.2<sup>-</sup>(-1)</td></tr><tr><td>OH</td><td><sup>2</sup>Π<sub>3/2</sub> J = 5/2 ← 3/2</td><td>F = 2<sup>-</sup> → 2<sup>+</sup></td><td>2514.2987</td><td>+2.1</td><td>1.4<sup>-</sup>(-2)</td><td rowspan="3">121</td><td rowspan="3">SOFIA/GREAT</td><td rowspan="3">0.70</td><td rowspan="3">12′′</td><td rowspan="3">(2)</td></tr><tr><td rowspan="3">HDO</td><td rowspan="3">1<sub>1,1</sub>-0<sub>0,0</sub></td><td>F = 3<sup>-</sup> → 2<sup>+</sup></td><td>2514.3167</td><td>0.0</td><td>1.4<sup>-</sup>(-1)</td></tr><tr><td>F = 2<sup>-</sup> → 1<sup>+</sup></td><td>2514.3532</td><td>-4.4</td><td>1.2<sup>-</sup>(-1)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">p-H<sub>2</sub><sup>18</sup>O</td><td rowspan="2">1<sub>1,1</sub>-0<sub>0,0</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Notes.** Frequencies, Einstein coefficients, and upper-level energies are from the JPL (Pickett et al. 1998) database. The <sup>18</sup>OH line is only available for G34.26+0.15. HDO and the H<sub>2</sub><sup>18</sup>O line has only been used for G34.26+0.15. Molecules below the horizontal line are archival data. <sup>(*)</sup>Velocity offset of the hfs components. <sup>(†)</sup>The Einstein coefficients are listed in the form of A(B) = A × 10<sup>B</sup>. <sup>(‡)</sup>This line has been observed in Cycle 4 towards the source G35.20-0.7 with a main beam efficiency of 0.70.

**References.** (1) This work, (2) Wiesemeyer et al. (2016), (3) Coutens et al. (2014), (4) Liu (2017), (5) Flagey et al. (2013).

Herschel Science Archive and processed the data with the IRAM/GILDAS software. We averaged the horizontal and vertical polarisations and summed up the spectra from the three different local oscillator (LO) tunings that have been used to identify line contamination from the image band. We fitted the spectral baseline over a narrow velocity range around the rest velocities (v<sub>lsr</sub>) of the sources and assumed a sideband gain ratio of unity to correct for the double-sideband receiver. We then used a beam efficiency of 74% (Roelfsema et al. 2012) to convert the spectra to a main beam brightness temperature scale.

For G327.29-0.6, we used the HDO (J<sub>K,A</sub> = 1<sub>1,1</sub>-0<sub>0,0</sub>) line at 893.6 GHz from the WISH (van Dishoeck et al. 2021) programme also observed with the Herschel/HIFI instrument. The data were processed in the same way as above, with the exception that in this case a single frequency setup was used to cover the line.

## 2.4. Archival APEX data

Towards the source G351.58-0.4, the complementary data for the HDO (J<sub>K,A</sub> = 1<sub>1,1</sub>-0<sub>0,0</sub>) line at 893.6 GHz were observed with the Atacama Pathfinder Experiment Telescope (APEX, Güsten et al. 2006) using the CHAMP+ receiver (Liu 2017).

# 3. Results

## 3.1. The OD $^2\Pi_{3/2}$ ($J = 5/2-3/2$) transition

### 3.1.1. Detection rates

As an example, in Fig. 2 we show the line-to-continuum ratio spectrum (T<sub>1</sub>/T<sub>c</sub>) of the OD $^2\Pi_{3/2}$ ground-state (J = 5/2-3/2) transition at 1391.5 GHz for the quiescent source G23.21-0.3.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Example of the OD $^2\Pi_{3/2}$ ground-state (J = 5/2-3/2) absorption towards the quiescent source G23.21-0.3. The spectra are normalised to the (single side-band) continuum level, and thus show the line to continuum ratio. The rest frequency is set to the strongest hyperfine transition, F = 7/2<sup>-</sup>5/2<sup>+</sup>, at 1391494.70 MHz. The blue dotted line marks the systemic LSR velocity of the source, and the grey dotted line shows the value of 1. The blue line shows a single component Gaussian fit to the spectrum. The positions and the relative strength of the hyperfine components are shown in green.</figcaption></figure>

Spectra for the rest of our sources are shown in Fig. 3. The lines typically exhibit a Gaussian profile that is a blend of the six hfs components of the Λ-doublet transition spanning a velocity range of ~4.2 km s<sup>-1</sup> (see also Table 2).

The OD line is detected towards nine of the 13 targets (69%) on the 3.5-22σ level. In many cases, it appears as a strong and narrow absorption feature. We first fitted the spectra with a
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Overview of the OD $^2\Pi_{3/2}$ ground-state (J = 5/2–3/2) transition towards the sample. The spectra are normalised to the continuum level and thus show the line-to-continuum ratio. The rest frequency is set to the strongest hyperfine transition, F=7/2⁻–5/2⁺, at 1391494.70 MHz. The vertical blue dotted line marks the systemic LSR velocity of the source (Table 1), horizontal grey dotted line shows the intensity ratio of 1. The light blue line shows a single or a two-component Gaussian fit to the spectrum. The positions and the relative strength of the hyperfine components are shown in green. The red spectra correspond to the HDO (J$_{K,A}$ = 1$_{1,1}$–0$_{0,0}$) line at 893.6 GHz from Herschel/HIFI, the orange spectrum shows the same line observed with the APEX telescope. The source name is labelled in the lower right corner of each panel.

single or a two-component Gaussian to extract the line parame-
ters and evaluate the detection rate. The fitted line parameters
are summarised in Table 3, which also includes the measured
continuum emission near the OD line frequency accounting for
the double-sideband reception of the signal. Among the non-
detections, some sources may show a weak feature close to our
signal-to-noise criterion; for example, towards W33A an absorp-
tion feature around 2σ<sub>rms</sub> is observed close to the expected LSR
velocity (v<sub>lsr</sub>). In addition, towards W28A the features around
the expected v<sub>lsr</sub> could also be consistent with a weak and broad
absorption component and also a line in emission at the 2σ<sub>rms</sub>
level. The line profiles of these components would resemble
those of the <sup>2</sup>Π<sub>3/2</sub> OH line towards this source analysed in
Leurini et al. (2015), although we caution that in our data the
baseline quality for this source is inferior to that of the other
sources.

3.1.2. The velocity components of OD absorption

In Figs. 2 and 3, we show the $v_{\text{lsr}}$ of the sources from Table 1. The OD $^2\Pi_{3/2}$ ($J = 5/2-3/2$) absorption is typically redshifted compared to the source's $v_{\text{lsr}}$, which is clearly the case for G34.26+0.15, G351.58-0.4, G10.47+0.03 and W31C, corresponding to 44% of the detections. For these sources, the HDO ($J_{\text{K,A}}=1_{1,1}-0_{0,0}$) line is also seen similarly redshifted, as well as the NH$_3$ ($J = 3_{2+}-2_{2-}$) transition from Wyrowski et al. (2016). The rest velocity at the peak of the absorption is on average ~1.5 km s$^{-1}$ redshifted for both the OD $^2\Pi_{3/2}$ ($J = 5/2-3/2$) and NH$_3$ ($J = 3_{2+}-2_{2-}$) lines.

The OD $^2\Pi_{3/2}$ ($J = 5/2-3/2$) absorption feature coincides well with the source systemic velocity for G35.20-0.7, W51e2 and G327.29-0.6, while the NH$_3$ ($J = 3_{2+}-2_{2-}$) line is detected as redshifted by Wyrowski et al. (2016). Considering the errors of the fit, which are larger by several factors for the OD line

compared to the NH₃ line, we find that the v<sub>lsr</sub> of OD and NH₃
are also consistent for these sources.

The OD $^2\Pi_{3/2}$ ($J = 5/2-3/2$) absorption is seen blueshifted
as compared to the source's $v_{\text{lsr}}$ towards W49N, often interpreted
as a P-Cygni profile due to expanding gas. The same kinematic
feature has already been noted by Wyrowski et al. (2016) using
the NH$_3$ ($J = 3_{2+}-2_{2-}$) line, and it is likely related to the expand-
ing shell structure of this source (Peng et al. 2010). The OD
$^2\Pi_{3/2}$ ($J = 5/2-3/2$) absorption seen towards W49N needs to be
fitted with two velocity components; one corresponds to the typi-
cally observed narrow absorption feature at the systemic velocity
like for the other sources, while the other absorption feature is
seen over a broad velocity range of $\Delta v = 11 \text{ km s}^{-1}$. This could
be explained if it was associated with outflowing gas at lower
velocities than the source's $v_{\text{lsr}}$.

We detect and resolve two velocity components in absorp-
tion towards W51e2, one at the $v_{\text{lsr}}$ of the main cloud, and a
second component associated with the more diffuse foreground
cloud unrelated to W51 Main at ~70 km s⁻¹ (Mookerjea et al.
2014). This is surprising because this gas has been considered to
originate from a layer of diffuse gas in front of the W51 Main
complex, unrelated to its star formation activity and should have
a much lower column density.

It is clear that OD is predominantly found towards massive
envelopes of star forming regions in various evolutionary stages.
However, considering its detection both in the (presumably)
more diffuse medium of the 70 km s⁻¹ component of W51e2 and
in outflowing gas from W49N, it appears that it could also be
present in sources that have different physical conditions.

3.1.3. The velocity structure of the line

Due to the hyperfine structure of the line, the Gaussian fit
overestimates the velocity dispersion. Therefore, we also fit the
<table><caption>**Table 3.** Fit results and gas parameters estimated from the $^2\Pi_{3/2}$ OD $J=5/2-3/2$ line.</caption><thead><tr><th>Source</th><th>T<sub>c</sub><sup>OD (a)</sup><br>(K)</th><th>σ<sub>rms</sub><br>(K)</th><th>Method<sup>(b)</sup></th><th>T<sub>mb</sub><sup>(c)</sup><br>(K)</th><th>v<sub>lsr</sub><sup>(d)</sup><br>(km s<sup>-1</sup>)</th><th>Δv<sup>(c)</sup><sub>FWHM</sub><br>(km s<sup>-1</sup>)</th><th>T<sub>1</sub>/T<sub>c</sub><sup>(e)</sup><br>(K)</th><th>τ<sub>tot</sub><sup>(e)</sup></th><th>N<sub>tot</sub><sup>(e)</sup><br>×10<sup>12</sup> (cm<sup>-2</sup>)</th></tr></thead><tbody><tr><td>G10.47+0.03</td><td>3.7</td><td>0.17</td><td>Gaussian<br>hfs</td><td>-0.27 ± 0.06</td><td>68.13 ± 0.37<br>68.90 ± 0.37</td><td>4.45 ± 0.69<br>3.63 ± 0.85</td><td>0.07 ± 0.05</td><td>0.08 ± 0.04<br>0.08</td><td>5.30 ± 2.65</td></tr><tr><td>G23.21-0.3</td><td>0.8</td><td>0.06</td><td>Gaussian<br>hfs</td><td>-0.47 ± 0.05</td><td>77.01 ± 0.13<br>77.80 ± 0.13</td><td>3.49 ± 0.27<br>2.20 ± 0.30</td><td>0.59 ± 0.09</td><td>0.89 ± 0.05<br>0.89</td><td>35.76 ± 2.01</td></tr><tr><td>G34.26+0.15</td><td>6.5</td><td>0.09</td><td>Gaussian<br>hfs</td><td>-1.31 ± 0.05</td><td>60.31 ± 0.05<br>61.10 ± 0.05</td><td>4.18 ± 0.12<br>3.24 ± 0.14</td><td>0.20 ± 0.01</td><td>0.23 ± 0.01<br>0.23</td><td>13.61 ± 0.59</td></tr><tr><td>G34.41+0.2</td><td>1.4</td><td>0.10</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>G35.20-0.7</td><td>1.7</td><td>0.07</td><td>Gaussian<br>hfs</td><td>-0.60 ± 0.06</td><td>33.13 ± 0.11<br>33.90 ± 0.12</td><td>3.70 ± 0.29<br>2.60 ± 0.36</td><td>0.35 ± 0.04</td><td>0.44 ± 0.03<br>0.44</td><td>20.89 ± 1.42</td></tr><tr><td>W28A</td><td>3.9</td><td>0.11</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>W31C</td><td>6.5</td><td>0.07</td><td>Gaussian<br>hfs</td><td>0.19 ± 0.06<br>-0.61 ± 0.16</td><td>-8.29 ± 0.51<br>-0.39 ± 0.17<br>-2.65 ± 0.66<br>-0.30 ± 0.66</td><td>4.96 ± 1.23<br>3.83<sup>(†)</sup><br>9.10 ± 2.19<br>6.62 ± 2.19</td><td>0.10 ± 0.01</td><td>0.10 ± 0.01</td><td>12.09 ± 1.21</td></tr><tr><td>W33A</td><td>1.7</td><td>0.06</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.1</td></tr><tr><td>W49N</td><td>9.5</td><td>0.10</td><td>Gaussian<br>hfs</td><td>-1.13 ± 0.09<br>-0.38 ± 0.04</td><td>10.78 ± 0.17<br>3.30 ± 0.83<br>3.97 ± 1.20<br>11.50 ± 0.18</td><td>5.79 ± 0.34<br>12.26 ± 1.06<br>11.70 ± 1.52<br>5.28 ± 0.45</td><td>0.12 ± 0.01<br>0.04 ± 0.01</td><td>0.13 ± 0.01<br>0.04 ± 0.01<br>0.13<br>0.04</td><td>27.78 ± 2.14<br>3.86 ± 0.96</td></tr><tr><td>W51e2</td><td>7.6</td><td>0.13</td><td>Gaussian<br>hfs</td><td>-1.06 ± 0.13<br>-0.66 ± 0.09</td><td>68.06 ± 0.14<br>57.40 ± 0.31<br>68.80 ± 0.11<br>58.2 ± 0.31</td><td>3.63 ± 0.30<br>7.10 ± 0.74<br>1.63 ± 0.32<br>6.74 ± 0.81</td><td>0.14 ± 0.02<br>0.09 ± 0.02</td><td>0.15 ± 0.02<br>0.09 ± 0.02<br>0.15<br>0.09</td><td>4.47 ± 0.60<br>11.08 ± 2.46</td></tr><tr><td>W51d</td><td>4.7</td><td>0.08</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>G327.29-0.6</td><td>5.6</td><td>0.08</td><td>Gaussian<br>hfs</td><td>-1.32 ± 0.08</td><td>-44.40 ± 0.09<br>-43.60 ± 0.9</td><td>5.22 ± 0.23<br>4.43 ± 0.25</td><td>0.24 ± 0.01</td><td>0.27 ± 0.01<br>0.27</td><td>21.85 ± 0.81</td></tr><tr><td>G351.58-0.4</td><td>1.7</td><td>0.11</td><td>Gaussian<br>hfs</td><td>-1.23 ± 0.09</td><td>-94.69 ± 0.11<br>-93.90 ± 0.11</td><td>4.36 ± 0.25<br>2.94 ± 0.27</td><td>0.72 ± 0.08</td><td>1.29 ± 0.05<br>1.29</td><td>69.27 ± 2.68</td></tr></tbody></table>

**Notes.** (a) Main-beam brightness temperature of the single-sideband continuum at 1.391 THz (Rayleigh-Jeans scale). (b) The fit corresponds to a one- or two-component Gaussian fit to the spectra, while the hfs method takes into account the relative strengths and velocity pattern of the hyperfine structure lines. (c) The errors are propagated from the HFS fit errors from CLASS using the respective formulae. (d) The errors are constrained by the fitting procedure in CLASS, assuming a normal distribution of uncertainties. (e) The errors are propagated from the fit errors and the noise of the data ($\sigma_{rms}$) using standard formulation. (f) The line width is fixed to that of the HDO line shown in Fig. 3.

spectra using the GILDAS/CLASS HFS fitting tool that takes into account the hfs splitting and is based on the assumption that the hfs components have the same excitation temperature. This fit constrains the $v_{lsr}$ and the line width ($\Delta v$) relatively well; however, the line opacity ($\tau$) is poorly constrained due to the blending of the components causing large errors in the opacity estimate. To overcome this, we used the line opacity estimated from the line-to-continuum ratio discussed in Sect. 3.1.4 and fixed this value for the HFS fit. The values obtained this way are reported in Table 3.

Comparing the results of the HFS fit to those of a simple Gaussian fit, the HFS fitting reveals that the blending of the six hfs components broadens the line widths considerably, typically by up to a factor of two compared to a simple Gaussian fit to the blended features. In general, the OD absorption is a typically narrow feature with a line width, $\Delta v$, between 2.2 and 11.1 km s$^{-1}$, a mean of 4.2, and a median of 3.2 km s$^{-1}$. The broadest component is detected towards W49N, and it is likely associated with the outflowing material that is also seen in the $^2\Pi_{3/2}$ excited OH line (Sect. 3.2). The emission component seen towards W31C is somewhat broader or blueshifted than the absorption component, suggesting either a higher level of turbulent broadening compared to the absorption feature or a different velocity component.

This blueshifted component at approximately $-8$ km s$^{-1}$ has also been seen in the $\mathrm{o-NH_3}$ ($J = 2_0-1_0$) line at 1214.9 GHz (Persson et al. 2010).

The line width from the OD absorption line is, on average, $4.9$ km s$^{-1}$ with a median of $3.6$ km s$^{-1}$. This is considerably narrower than that of the $\mathrm{C^{17}O}$ ($J = 3-2$) line from Giannetti et al. (2014), which is found to be, on average, $5.7$ km s$^{-1}$ (with a median of $5.8$ km s$^{-1}$). The dispersion on the line-width estimates from the $\mathrm{C^{17}O}$ ($J = 3-2$) line exhibit, however, a considerably smaller dispersion of $1.7$ km s$^{-1}$, while the OD line widths have a standard deviation of $3$ km s$^{-1}$ over the sample. This suggests that the OD $^2\Pi_{3/2}$ ($J = 5/2-3/2$) line tends to trace a somewhat more quiescent gas compared to that of the $\mathrm{C^{17}O}$ ($J = 3-2$) line.

### 3.1.4. OD column density

The OD $^2\Pi_{3/2}$ ($J = 5/2-3/2$) line is detected in absorption towards all sources, and an emission component is only seen towards W31C with a $3.5 \sigma_{rms}$ detection measured over the peak brightness temperature of the line. In contrast to the low-mass protostar, IRAS 16293-2422 (Parise et al. 2012), the OD absorption towards massive Galactic star forming regions does not
<table><caption>**Table 4.** Fit results of the $^2\Pi_{1/2}$ OH $J = 3/2-1/2$ line.</caption><thead><tr><th>Source</th><th><i>T</i><sub>c</sub><sup>OH (a)</sup><br>(K)</th><th><i>σ</i><sub>rms</sub><br>(K)</th><th><i>T</i><sub>mb</sub><sup>(b)</sup><br>(K)</th><th><i>v</i><sub>lsr</sub><sup>(b)</sup><br>(km s<sup>-1</sup>)</th><th>Δ<i>v</i><sup>(b)</sup><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>G10.47+0.03</td><td>5.49</td><td>0.06</td><td>1.2 ± 0.2<br>-0.8 ± 0.1</td><td>64.3 ± 0.4<br>49.5 ± 1.7</td><td>4.3 ± 1.1<br>26.9 ± 3.0</td></tr><tr><td>G34.26+0.15<sup>(†)</sup></td><td>8.66</td><td>0.14</td><td>1.6 ± 0.2<br>-1.1 ± 0.3</td><td>58.1 ± 0.3<br>62.5 ± 0.3</td><td>4.2 ± 0.7<br>1.5 ± 0.6</td></tr><tr><td>G35.20-0.7</td><td>2.02</td><td>0.15</td><td>2.1 ± 0.2</td><td>34.5 ± 0.4</td><td>9.7 ± 0.9</td></tr><tr><td>W28A</td><td>6.80</td><td>0.11</td><td>5.7 ± 2.1<br>4.0 ± 1.5</td><td>26.0 ± 6.2<br>12.4 ± 1.3</td><td>34.2 ± 8.5<br>13.7 ± 5.8</td></tr><tr><td>W31C</td><td>8.79</td><td>0.12</td><td>3.5 ± 0.3<br>2.0 ± 0.3</td><td>-3.3 ± 0.5<br>4.2 ± 1.2</td><td>7.8 ± 1.0<br>14.9 ± 2.9</td></tr><tr><td>W33A</td><td>2.32</td><td>0.13</td><td>1.5 ± 0.1</td><td>40.7 ± 0.4</td><td>16.3 ± 1.2</td></tr><tr><td>W49N</td><td>12.35</td><td>0.07</td><td>4.8 ± 0.4<br>3.8 ± 0.3<br>2.9 ± 0.2</td><td>19.8 ± 0.4<br>6.7 ± 1.2<br>13.9 ± 3.8</td><td>10.2 ± 1.1<br>12.8 ± 1.2<br>68.6 ± 3.8</td></tr><tr><td>G327.29-0.6</td><td>5.53</td><td>0.06</td><td>—</td><td>—</td><td>—</td></tr><tr><td>G351.58-0.4</td><td>6.17</td><td>0.10</td><td>0.7 ± 0.2</td><td>-86.0 ± 1.9</td><td>15.9 ± 3.9</td></tr></tbody></table>

**Notes.** $^{(a)}$Main-beam brightness temperature of the single-sideband continuum at 1.838 THz (Rayleigh-Jeans scale). The velocity resolution is 0.75 km s$^{-1}$. $^{(b)}$The errors are constrained by the fitting procedure in CLASS. $^{(†)}$The hfs fit is for the 1834.8 GHz component of the Λ-doublet from the image band.

appear to be highly saturated; instead, it is typically optically thin. The deepest absorption with $T_1/T_c$ ratios of 0.6 and 0.7 is seen towards G23.21-0.3 and G351.58-0.4, respectively.

The deepest absorption together with the smallest line-widths are observed towards sources that are typically the youngest, like G23.21-0.3, or that are still dominated by a massive, but relatively cold envelope, such as G351.58-0.4, G327.29-0.6. The total optical depth is $\tau_{tot} = \sum_i \tau_i$, where $i$ corresponds to the hfs component of the line, and we estimated it using the formula $\tau_{tot} = -\ln(1 - \frac{T_1}{T_c})$, where $T_1$ corresponds to the main beam brightness temperature of the line and $T_c$ is the continuum temperature accounting for the double-sideband reception of the signal. We find that the values for $\tau_{tot}$ range between 0.04 and 1.29, with a median of 0.15 and a mean of 0.34, implying that the absorption is typically optically thin. The largest optical depths ($\tau_{tot}$) measured towards G351.58-0.4 and G23.21-0.3 are 1.29 and 0.89, respectively. The smallest values are seen towards the broad absorption component of W49N, and the 70 km s$^{-1}$ diffuse cloud component of W51e2.

Adopting the collisional rate coefficients of OH also for OD (see Sect. 4), we find that the $^2\Pi_{3/2}$ OD $J = 5/2-3/2$ line has a very high critical density, $n_{crit}$ of $> 10^8$ cm$^{-3}$. Such a high volume density is unlikely to be present in the gas in which the OD absorption originates. Therefore, LTE conditions are unlikely to apply. Moreover, the upper level energy ($h\nu/k = 67$ K) largely exceeds the temperature of the cosmic microwave background. Therefore, it is safe to assume that for all sources the vast majority of OD molecules are in the rotational ground state. The total column density of the molecule can then be estimated as (e.g. Comito et al. 2003) $N_{tot} = \sum_i \frac{8\pi v_i^3}{A_{ul,i}c^3} \Delta v_i \frac{g_{l,i}}{g_{u,i}} \tau_i$, where $v_i$ is the transition frequency and $A_{ul,i}$ the Einstein-coefficient of transition $i$. $\Delta v_i$ is the equivalent width of the line, which relates to its full width at half maximum $\Delta v_{FWHM}$ reported in Table 4 via $\Delta v = \Delta v_{FWHM} \sqrt{\pi/(4 \ln 2)}$. $g_{l,i}$ and $g_{u,i}$ are the lower and upper state degeneracies, and $\tau_i$ is the optical depth at the line centre of the $i$ hfs component. The summation extends over the observed Λ-doublet component and applies a correction for the OD population in the other one so as to obtain $N_{tot}$. We

find values between $3.9 \times 10^{12}$ and $6.9 \times 10^{13}$ cm$^{-2}$, where the lowest values are typically measured towards the more evolved sources hosting UCHII regions, and the highest values towards the sources that are dominated by a significant amount of cold dust corresponding to quiescent regions and hot cores. This conclusion is fairly robust against our assumption of a complete ground-state population: we repeated the analysis assuming an excitation temperature of 20 K instead, and the deviation in the deduced $N_{tot}$ is within or close to our fit errors.

### 3.2. The $^2\Pi_{1/2}$ OH ($J = 3/2-1/2$) transition

We observed the rotationally excited $^2\Pi_{1/2}$ OH ($J = 3/2-1/2$) transition at 1834.7 GHz towards ten sources of the sample, and show both Λ-doublet transitions in Fig. 4. Two sources without an OD detection, W51d and G31.41+0.2, as well as the quiescent source, G23.21-0.3, and the UCHII region, W51e2, have not been observed in this line. Furthermore, the excited OH line is not detected towards G351.58-0.4 and G327.29-0.6. Towards W33A, only one of the Λ doublets falls in the band.

The remaining sources show relatively bright emission in the $^2\Pi_{1/2}$ OH ($J = 3/2-1/2$) line. The Λ doublets have three hfs components (Table 2), the 1834.7 GHz line has an isolated hfs component separated by 11.5 km s$^{-1}$, which appears as a distinct feature in the spectra. The largest velocity shift for the other Λ-doublet component at 1837.8 GHz is separated by 1.9 km s$^{-1}$, so all its three hfs components appear blended in the spectra.

The OH line profiles frequently exhibit a broad velocity component likely associated with outflowing high-velocity gas, which is very pronounced in W49N, W28A, and W33A. Towards G10.47+0.03, the outflow component appears to be in absorption. Previously, outflowing gas has been found in excited OH towards the high-mass star forming region W3 IRS 5, where Wampfler et al. (2011) also detected both the envelope and outflow components. Our previous SOFIA/GREAT observations of this excited OH line did not have a high enough signal-to-noise ratio (S/N) to be sensitive to weaker broad components (Csengeri et al. 2012). Low spectral resolution observations with Herschel similarly suggest that the exited OH line can be
<table><caption>**Table 5.** Fit parameters from RATRAN modelling of the OH and OD lines.</caption><thead><tr><th>Source</th><th><i>d</i><br>(kpc)</th><th><i>L</i><sub>bol</sub><br>(L<sub>☉</sub>)</th><th><i>R</i><sub>out</sub><br>(pc)</th><th>α<sub>n</sub></th><th><i>n</i><sub>1 pc</sub><br>(10<sup>3</sup> cm<sup>-3</sup>)</th><th>Δ<i>v</i><br>(km s<sup>-1</sup>)</th><th><i>X</i>(OH)</th><th><i>X</i>(OD)</th><th><i>D</i><sup>OD</sup><sub>out</sub></th><th><i>D</i><sup>HDO</sup><sub>out</sub></th><th>Ref.<br>for <i>D</i><sup>HDO</sup></th></tr></thead><tbody><tr><td>G23.21–0.3</td><td>4.6</td><td>1.3 × 10<sup>4</sup></td><td>1.8</td><td>−2.0</td><td>4.5</td><td>1.0</td><td>—</td><td>8.0 × 10<sup>-10</sup></td><td>—</td><td>—</td><td>—</td></tr><tr><td>G34.26+0.2</td><td>1.6</td><td>4.6 × 10<sup>4</sup></td><td>0.8</td><td>−1.8</td><td>10</td><td>2.4</td><td>1.3 × 10<sup>-8</sup></td><td>6.2 × 10<sup>-11</sup></td><td>4.8 × 10<sup>-3</sup></td><td>5.1 × 10<sup>-3</sup><br>1.0–3.1 × 10<sup>-3</sup><br>1.9–4.9 × 10<sup>-4</sup></td><td>(a)<br>(b,c)<br>(d)</td></tr><tr><td>G34.41+0.2</td><td>5.4</td><td>6.7 × 10<sup>4</sup></td><td>2.5</td><td>−2.0</td><td>7</td><td>3.2</td><td>—</td><td>&lt;8.0 × 10<sup>-11</sup></td><td>—</td><td>—</td><td>—</td></tr><tr><td>G35.20–0.7</td><td>2.2</td><td>2.5 × 10<sup>4</sup></td><td>1.5</td><td>−1.6</td><td>5.5</td><td>1.5</td><td>&gt;2.0 × 10<sup>-8</sup></td><td>3.0 × 10<sup>-10</sup></td><td>&lt;1.5 × 10<sup>-2</sup></td><td>—</td><td>—</td></tr><tr><td>G327.29–0.6</td><td>3.1</td><td>8.2 × 10<sup>4</sup></td><td>2.0</td><td>−1.9</td><td>10</td><td>2.3</td><td>&lt;5 × 10<sup>-9</sup></td><td>9.75 × 10<sup>-11</sup></td><td>&gt;2.0 × 10<sup>-2</sup></td><td>5.9–6.4 × 10<sup>-3</sup></td><td>(e)</td></tr><tr><td>G351.58–0.4</td><td>6.8</td><td>2.4 × 10<sup>5</sup></td><td>1.8</td><td>−1.9</td><td>15</td><td>1.5</td><td>7.8 × 10<sup>-9</sup></td><td>1.95 × 10<sup>-10</sup></td><td>2.5 × 10<sup>-2</sup></td><td>0.6–1.4 × 10<sup>-3</sup></td><td>(e)</td></tr></tbody></table>

**Notes.** Outer radius, $R_{\text{out}}$, density power-law index, $\alpha_n$, and the density at 1 pc, $n_{1\text{pc}}$, were constrained by the 870 µm ATLASGAL dust continuum emission and are adopted from Wyrowski et al. (2016), as are the distances.

**References.** (a) This work. (b) Coutens et al. (2014). (c) Kulczak-Jastrzębska (2017). (d) Liu et al. 2013. (e) Liu (2017).

frequently associated with emission from the high velocity outflowing gas (cf. Yang et al. 2018). The high-velocity outflowing gas of W28A has been studied in detail, including the OH data shown here, in Leurini et al. (2015).

In addition to the broad velocity component, a narrow component associated with the envelope is always present. Furthermore, the line profiles towards G34.26+0.15 suggest that both the emission and absorption components are present in the excited OH line, similarly to what is seen in other high density gas tracers (cf. Wyrowski et al. 2016). The spectra towards G10.47+0.03 are a mixture of a narrow emission component and a broader absorption component from the outflowing gas, as discussed above. The double-peaked profile towards W49N suggests absorption from lower-excitation gas in the envelope.

Similarly to our previous work in Csengeri et al. (2012), we used the HFS method of CLASS to fit the lines. We did this separately for the Λ doublets in the signal and image side band, respectively, and list the average of these results in Table 4 and display them in Fig. 4. Except for W49N, where three components were needed, typically one narrow and a broader component was sufficient to fit the data. The brightest emission is seen towards sources that appear typically weak in the ground-state OD line.

The average FWHM line width for the $^2\Pi_{1/2}$ OH ($J=3/2-1/2$) line is 10.3 km s$^{-1}$ for all the sample considering the narrowest components, where two velocity components were fitted. This is more than two times larger than that of the OD $^2\Pi_{3/2}$ ($J=5/2-3/2$) line, where for the sources also observed in the $^2\Pi_{1/2}$ OH ($J=3/2-1/2$) line we find an average line-width of 3.6 km s$^{-1}$.

# 4. Line modelling with RATRAN

## 4.1. The physical structure of the envelope

To estimate the abundance distribution of OD and to compare it to that of OH throughout the envelope, we use the non-LTE radiative transfer code RATRAN (Hogerheijde & van der Tak 2000). Following the work of Wyrowski et al. (2016), we constrained the physical structure of the envelope (such as size and the density profile) with the 870 µm dust emission as observed in the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL, Schuller et al. 2009). This implies that we assume that the densities are high enough for the gas and dust to be in thermal equilibrium, which is a reasonable assumption for high-mass

protostellar envelopes embedded in massive clumps with high densities. The corresponding parameters of the physical properties are listed in Table 5 for sources with a detection in the ground-state OD line. This allowed us to model six sources.

Our models use a spherical geometry with a power-law profile for the density, $n$, as a function of radius, $r$, given by $n(r) \sim r^{\alpha_n}$, and a temperature, $T$, profile of $T(r) \sim r^{-0.4}$. The line width is obtained from the $^{17}$O (3–2) data for these sources (Giannetti et al. 2014). The velocity gradient through the envelope is constrained by $\text{NH}_3$ observations and the models of Wyrowski et al. (2016). The physical structure of the model for G34.26+0.15 is illustrated in Fig. 5, where the bottom panel shows the modelled abundance profiles for OH and OD (see Sect. 4.2). We used the recent collisional rate coefficients from Kłos et al. (2017) for OH for collisions with $\text{H}_2$, and adopted these coefficients for OD (see Appendix D for the validity of this approach). We then varied the OH and OD abundances ($X(\text{OH})$ and $X(\text{OD})$), respectively.

Radiative heating from (high-mass) protostars leads to a radial increase in temperature towards the centre of the envelope (e.g. Wolfire & Cassinelli 1986; Goldreich & Kwan 1974). At low temperatures ($\sim\sim 20$ K), the chemistry is determined by grain surface processes. Triggered by the ignited young stellar object, molecules are desorbed. This greatly changes gas-phase abundances, which are further modified by gas-phase reactions. More importantly, deuterium fractionation of several species is expected to be different in gas-phase reactions compared to grain surface processes (e.g. Roberts & Millar 2007), leading to different molecular abundances in the inner (hot) and outer (cold) regions of envelopes. Since the grain mantles are expected to be mostly composed of water ice, the typical temperature (and thus radius) determining the change of abundances is related to the water evaporation temperature around 100 K (Fraser et al. 2001).

Using 1D, non-LTE radiative transfer modelling, our aim is to estimate the OD abundance in the envelope and to probe what fraction of the OD originates from the hot inner regions and the cold outer envelope. Then, assuming that OH and OD are part of the same gas, we constrained their abundances throughout the envelope and thereby determined the deuterium fractionation by $D^{\text{OD}} = X(\text{OD})/X(\text{OH})$.

For our models, the ground-state $^2\Pi_{3/2}$ OH line seen in absorption should, in principle, provide the best constraint on the OH abundance in the outer envelope. However, as discussed in detail by Wiesemeyer et al. (2016), this line is largely optically thick and is severely blended with absorption from dif-fuse gas along the line of sight (cf. Wiesemeyer et al. 2012, 2016). Therefore, we can only check whether our models are consistent
<table><caption>**Table 6.** Summary of RATRAN model parameters and results towards G34.26+0.15.</caption><thead><tr><th>Mod.</th><th>X(OD)<br>profile</th><th>r<sub>j</sub><br>(")</th><th>X(OH)<br>model range (result)</th><th>X(OD)<sup>(†)</sup><br>model range (result)</th><th>D<sub>in</sub><sup>OD</sup></th><th>D<sub>out</sub><sup>OD</sup></th></tr></thead><tbody><tr><td>A</td><td>Constant</td><td>—</td><td>10<sup>-10</sup>–1 × 10<sup>-8</sup> (1.3 × 10<sup>-8</sup>)</td><td>1.0 × 10<sup>-14</sup>–2.5 × 10<sup>-9</sup> (6.2 × 10<sup>-11</sup>)</td><td>—</td><td>10<sup>-4</sup>–5 × 10<sup>-1</sup> (4.8 × 10<sup>-3</sup>)</td></tr><tr><td>B</td><td>T<sub>j</sub> = 50 K</td><td>5.0"</td><td>1.3 × 10<sup>-8</sup></td><td>1.3 × 10<sup>-12</sup>–6.5 × 10<sup>-10</sup></td><td>10<sup>-4</sup>–5 × 10<sup>-1</sup> (≤ 5 × 10<sup>-3</sup>)</td><td>4.8 × 10<sup>-3</sup></td></tr><tr><td>C</td><td>T<sub>j</sub> = 100 K</td><td>1"0</td><td>1.3 × 10<sup>-8</sup></td><td>1.3 × 10<sup>-12</sup>–6.5 × 10<sup>-10</sup></td><td>10<sup>-3</sup>–5 × 10<sup>-1</sup>(‡)</td><td>4.8 × 10<sup>-3</sup></td></tr><tr><td>D</td><td>T<sub>j</sub> = 200 K</td><td>0"2</td><td>1.3 × 10<sup>-8</sup></td><td>1.3 × 10<sup>-12</sup>–6.5 × 10<sup>-10</sup></td><td>10<sup>-3</sup>–5 × 10<sup>-1</sup>(‡)</td><td>4.8 × 10<sup>-3</sup></td></tr></tbody></table>

**Notes.** Bold face values correspond to the best fit results. (†)X(OD) corresponds to X(OD)in for models B, C, and D. (‡)Due to the small source size, all these models are consistent with the data.

with these observations. However, towards G34.26+0.15 we obtained archival data from SOFIA covering the $^{2}\Pi_{3/2}$ $^{18}$OH ($J = 5/2-3/2$) ground-state transition at 2494.6950 GHz. This line shows an emission spike at ~50 km s$^{-1}$ due to a telluric line, which, against strong continuum sources, makes the calibration locally challenging (see Fig. 6 and Wiesemeyer et al. 2012 for details). Nevertheless, we also included the $^{18}$OH line in our models assuming an isotopic ratio of $^{16}$O/$^{18}$O = 400 (Wilson 1999).

The $^{2}\Pi_{1/2}$ OH ($J = 3/2-1/2$) excited-state transition has an upper energy level of 270 K compared to the 67 K of the $^{2}\Pi_{3/2}$ OD ($J = 5/2-3/2$) transition, and so it is clearly a more sensitive probe of the warmer inner regions. Nevertheless, we also consider this line to further constrain our models, allowing us to estimate the total OH column density in the envelope.

## 4.2. G34.26+0.15: Non-LTE modelling of OH, OD, HDO, and $H_2^{18}O$ using RATRAN

### 4.2.1. Modelling constraints

The source G34.26+0.15 exhibits a profile with redshifted self-absorption in the excited-state OH line, and its source structure is relatively well constrained. Hence, it is for this source that we can perform a more detailed modelling in order to test different OD abundance profiles. The simplest model (A) corresponds to a constant abundance profile through the envelope, while models B, C, and D introduce an abundance jump of OD above a temperature ($T_j$) of 50, 100, and 200 K, respectively (Fig. 5, bottom panel). OD originating from the warm inner envelope is expected to show up in emission in the blueshifted part of the spectrum, allowing us to constrain the deuterium fractionation in this source for both the inner, warm, and outer cold envelope (appearing in absorption). We list the model parameters, radius ($r_j$), and temperature ($T_j$) corresponding to the abundance jump, together with the range of abundances and the modelling results, in Table 6.

For model A, we created a grid of parameters for the OH abundance, X(OH), between $10^{-10}$ and $1 \times 10^{-8}$. We modelled the deuterium fractionation $D^{OD} = X(OD)/X(OH)$ between $10^{-4}$ and $5 \times 10^{-1}$. These ranges were sampled by 20 points on a logarithmic scale. For models B and C, we used the result from model A to fix the OH abundance and the deuterium fractionation in the outer envelope, and we varied $D^{OD}_{in}$ as listed in Table 6 in an inner radius between $r_j = 0''2$ and $5''$, which translates to physical scales of ~320 and 8000 au, respectively. After performing the calculations with RATRAN, we convolved the modelled spectra with the beam of the telescope at the respective frequencies shown in Table 2 and the hfs structure assuming LTE conditions, meaning we adopted the relative hfs intensities according to their rates of spontaneous emission or

absorption$^4$. The models have been compared to the observations and were evaluated by a reduced $\chi^2$ method.

For the other molecules shown in Fig. 6, we used the same physical model as for OH and OD, and we used the abundance of Wyrowski et al. (2016) for NH$_3$, and a constant abundance profile for HDO and $H_2^{18}O$. Since we used the revised distance estimate of 1.4 kpc based on maser parallax measurements (Kurayama et al. 2011) for G34.26+0.15, compared to Coutens et al. (2014), we find that we need to scale their abundances of HDO and $H_2^{18}O$ for the outer envelope by factors of 5 and 1.5, respectively, in order to obtain a good fit to the data$^5$.

### 4.2.2. Constant deuteration profile towards G34.26+0.15

Our first result for G34.26+0.15 is that already a constant OH and OD abundance profile throughout the entire envelope (model A) gives a good fit to all the lines shown in Fig. 6. Our best fit gives an OD abundance, X(OD), of $6.2 \times 10^{-11}$, and we find that an OH abundance, X(OH), of $1.3 \times 10^{-8}$ is consistent with the observed $^{2}\Pi_{3/2}$ OH ground-state absorption and the $^{2}\Pi_{1/2}$ excited-state line$^6$. This model corresponds to a constant deuterium fractionation of $4.8 \times 10^{-3}$ throughout the envelope. Our fit reproduces the 1834 GHz $\Lambda$-doublet pair of the $^{2}\Pi_{1/2}$ OH excited line relatively well, although it somewhat over estimates the absorption component. On the other hand, the fit to the 1837 GHz $\Lambda$ doublet is less good: it underestimates the emission component, and, as for the ground-state-line $\Lambda$-doublet pair, it overestimates the absorption component.

Using the same physical model, we also fitted the HDO and $p-H_2^{18}O$ ($J = 1_{1,1}-0_{0,0}$) lines with a constant abundance profile of X(HDO) = $1.6 \times 10^{-11}$ and X(H$_2^{18}O$) = $7.8 \times 10^{-12}$ (corresponding to X(H$_2$O) = $3.1 \times 10^{-9}$ assuming an isotopic ratio of $^{16}$O/$^{18}$O = 400). These lines were modelled in detail by Coutens et al. (2014), who adopted a distance of 3.3 kpc for the source, almost two times larger than what we used here based on Wyrowski et al. (2016). It is already interesting that we obtain a good fit to these lines without an abundance jump, suggesting that these observations are not sensitive to the hot inner regions. However, our models require an X(H$_2$O) more than three times lower compared to Coutens et al. (2014) for the outer envelope.

<footnote>$^4$ This approximation is not valid for the highly optically thick absorption seen in the ground-state $^{2}\Pi_{3/2}$ OH line.</footnote>

<footnote>$^5$ This gives a deuterium fractionation for water vapour in the outer envelope of $5.1 \times 10^{-3}$ instead of $1.5 \times 10^{-3}$ reported in Coutens et al. (2014).</footnote>

<footnote>$^6$ Our model over-predicts the 1.3 THz continuum for the OD $^{2}\Pi_{3/2}$ $J = 3/2-1/2$ transition by 0.65 K (corresponding to 10%), and it under-predicts the continuum by 0.9 K (10%) at 1.8 THz for the OH $^{2}\Pi_{1/2}$ $J = 5/2-3/2$ line. It shows a discrepancy of 3.0 K (corresponding to 50%) for the OH $^{2}\Pi_{3/2}$ $J = 3/2-1/2$ line at 2.5 THz.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Physical structure of models for G34.26+0.15. Bottom panel shows the OH (red) and OD (grey) abundance variations in the envelope for the different models. For a better visibility, the OD abundance in the inner envelope ($X(\text{OD})^{\text{in}}$) is shifted for the different models. All models test the same abundance range.</figcaption></figure>

In Fig. 6, we also show the observations and the model of the NH₃ *J* = 3<sub>2,1</sub>−2<sub>2,0</sub> line from Wyrowski et al. (2016) as a comparison. Our results therefore suggest a deuterium fractionation of *D*<sub>out</sub><sup>OD</sup> = 4.8 × 10<sup>−3</sup> for OD, very similar to that of HDO/H₂O of *D*<sub>out</sub><sup>HDO</sup> = 5.1 × 10<sup>−3</sup> constrained by our modelling here.

### 4.3. Impact of evaporation

Due to water ice evaporation, a change in the abundance profile and deuterium fractionation can be expected (see also Coutens

et al. 2014 and Sect. 5.3). Therefore, we also tested models introducing a change in deuterium fractionation at a given temperature (*T*<sub>J</sub>). In models *B*, *C*, and *D* we varied the temperature of the abundance jump (*T*<sub>J</sub> = 50, 100, 200 K) and varied the deuterium fractionation in the hot inner regions compared to the outer envelope (*D*<sub>in</sub><sup>OD</sup>) in a range of 10<sup>−4</sup>−5 × 10<sup>−1</sup>.

In Fig. 7, we show the results of our model *B* with *T*<sub>J</sub> = 50 K that introduces the abundance jump at a few thousands of au in the envelope, which is comparable to the physical scale of water evaporation found by Coutens et al. (2014). These models show a noticeable change in the blueshifted component of OD seen in emission as a result of the velocity gradient over this source. As we see below, given our sensitivity limit, this is the only model that can put a constraint on the OD column density in the inner regions at *r*<sub>in</sub> ~ 11 000 au. Our results suggest that the observations are consistent with a deuterium fractionation in the inner envelope (*D*<sub>in</sub><sup>OD</sup>) of up to ~5 × 10<sup>−3</sup>, which would imply that a jump of deuterium fractionation up to a factor of 2 could remain undetected in our observations. However, a strong deuterium enhancement towards the inner regions is clearly excluded by these models. This suggests that the OD abundance in the inner regions cannot be significantly higher than in the outer regions, and thus it does not favour an enhancement in deuterium fractionation in the inner regions compared to the cold envelope. A firmer constraint on the deuterium fractionation towards the innermost regions can only be obtained by observations of higher energy transitions of OD.

We notice that models *C* and *D* do not introduce any noticeable change in the model results compared to model *A*, and additional calculations show that for ‘jump’ temperatures *T*<sub>J</sub> ≥ 60 K, the hot emitting region becomes too small in angular size (≤3" corresponding to 4800 au) and their emission or absorption is heavily diluted in the SOFIA beam.

We also find that at abundances slightly higher than the values we investigate here, the emission of OH and OD becomes optically thick, and hence it cannot provide a view of the inner regions. Instead, as the higher energy transitions get excited, they should provide a better probe of the inner regions. This fact is somewhat compensated by the velocity structure of the source for G34.26+0.15 and G351.58−0.4 (Fig. 5), allowing us to probe deeper into the envelope due to the accelerated infall, which results in red-shifted velocities closer to the source. For G34.26+0.15 (in contrast to Coutens et al. (2014), whose best fitting models use a temperature of 200 K for an abundance jump), due to the revised distance for this source, our models at the same temperature correspond to a considerably more compact hotter region within 270 au. The Coutens et al. (2014) best fit corresponds to an extent of 1.7" (thus 6600 au) in their models, meaning to a value similar to that of our low *T*<sub>J</sub> models.

### 4.4. Other sources

Above, we discuss the radiative transfer modelling of G34.26+0.15 in detail. For the other sources, we refrained from such a comprehensive modelling in this work due to the lack of a similarly detailed knowledge about the source structure and the lack of sufficient ancillary data in the form of HDO and H<sub>2</sub><sup>18</sup>O lines. For some sources, though, we were able to rely on the source model of Wyrowski et al. (2016), and we used this as an input for RATRAN to obtain column density estimates primarily for OD and, where available, for OH based on the excited Λ doublets. All these models assume a constant OD abundance profile throughout the envelope similarly to model *A* for G34.26+0.15.
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Observed spectra towards G34.25+0.15 are shown in black, and the red line shows the fitted RATRAN model (model A). Dotted line marks the $v_{lsr}$ of the source. Magenta line shows the atmospheric transmission, where the absorption feature in the $^{18}$OH spectrum around 50 km s$^{-1}$ is due to a deliberately uncorrected telluric ozone feature.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Results of model B for $^2\Pi_{3/2}$ OD ($J = 5/2-3/2$) transition of G34.26+0.15 using $T_J = 50$ K. The different $D_{in}^{OD}$ values used for the models are shown in the figure legend.</figcaption></figure>

The parameters of Wyrowski et al. (2016) use the C$^{17}$O (3–2) measurements from Giannetti et al. (2014) to define the sources’ LSR velocity. They provide a good model for NH$_3$. We adopted the same line width for the OD ground-state absorption, which worked well for all sources, except for G351.58–0.4, for which we find a two times larger line width from the HFS fitting. Our models require this broader line component to provide a good fit to the OD spectrum (see also Appendix A). Similarly, G35.20–0.7 also requires a larger line width (of 2.6 km s$^{-1}$) than the

previous models; however, in this case we find that a constant OH abundance profile does not give a good fit to the excited OH line. The line profiles suggest self-absorption, and while our models reproduce the line profile, they typically underestimate the line intensity, suggesting that OH could be present in a higher abundance towards the warmer inner regions. We suggest that an increased OH column density in the inner envelope could provide a better fit to the observations. The reported OH column density is therefore an upper limit for the outer envelope, and by increasing its value even the excited OH line would appear only in absorption. Due to the limitations of the models for the OH abundance, we obtain the best constraints on the deuterium fractionation ($D^{OD}$) for G34.26+0.15. The results of these models are summarised in Table 5.

Altogether, the non-LTE radiative transfer modelling gives constraints on $X(\text{OD})$ towards six sources of the sample. These results agree within a factor of 3 with the ones calculated from the absorption features in Sect. 3.1. The largest differences are seen towards G23.21-0.3 and G327.29–0.6, where the modelling gives nearly three times higher OD abundances.

# 5. Discussion

## 5.1. OD absorption from the cold outer envelope

Absorption from the $^2\Pi_{3/2}$ OD ($J = 5/2-3/2$) ground-state line seems to be prevalent in high-mass star forming regions. It is, however, not exclusively associated with the high column density envelopes, but it is also observed towards various physical components within these regions. It is present in the broad absorption from the outflowing high-velocity gas towards W49N, as well as
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Total column density (N) of the ground-state $^2\Pi_{3/2}$ OD ($J = 5/2-3/2$) line versus the line width ($\Delta v$) obtained from the HFS fitting with a fixed total optical depth from the line-to-continuum ratio.</figcaption></figure>

in the more diffuse gas component seen, for example, towards the 70 km s⁻¹ cloud component around W51e2. In contrast to the low-mass protostar IRAS 16293–2422 (Parise et al. 2012), where the absorption is close to saturation, here the largest optical depths are seen towards the youngest, quiescent sources, and nearby hot cores, where a significant amount of cold dust is still present. These findings are consistent with the chemical models for OD by Croswell & Dalgarno (1985) that predict a high fractionation of OD in the low-density cold gas, while in the dense gas even higher fractionation could be expected depending on the local conditions, such as UV field and fractional ionisation.

Besides the potential outflow and diffuse gas origin, we focus here on the origin of OD within the high-mass protostellar envelopes, that is, the dense gas component. After accounting for the hfs splitting, the line widths we determine are narrow, which is consistent with an origin in the outer envelope where the gas is quiescent. Figure 8 shows a comparison between the measured line-width ($\Delta v$) from the OD absorption compared to the estimated total column density. The largest line widths are associated with a more turbulent component of the gas that seems to be less optically thick and have the smallest OD column density in the sample. The largest OD columns show narrow line widths, suggesting an origin from the dense gas of the envelope.

Comparing the $v_{lsr}$ and $\Delta v$ of the $^2\Pi_{3/2}$ OD ground-state line with those of the NH₃ ($J = 3_2-2_1$) absorption at 1.8 THz from Wyrowski et al. (2012, 2016), we notice that the line widths and the line shapes are in general similar (see Appendix B for more details). This further suggests a common origin in the cold absorbing layer of the envelope. Comparing with NH₃, in two cases we find a considerably larger line width in the OD line than in the NH₃ lines, suggesting a potential contribution from the more inner layers of the envelope. However, there are also other noticeable differences compared to the 1.8 THz NH₃ line. In contrast to the narrow OD absorption, NH₃ may exhibit velocity components that are not seen in the $^2\Pi_{3/2}$ OD ($J = 5/2-3/2$) line. For example, the broad absorption components seen in NH₃ towards the sources G327.29–0.6 and G351.58–0.4 are less pronounced in OD, suggesting either a chemically different environment, or different excitation conditions for this gas component.

The strongest confirmation that the $^2\Pi_{3/2}$ OD ($J = 5/2-3/2$) line originates from the cooler layers of the envelope is given by our detailed radiative transfer modelling of G34.26+0.15. We

find that the outer layers of the envelope at typically low temperatures ($T_d < 80–100$ K) and densities below $\sim 10^7$ cm⁻³ are responsible for the absorption feature (Fig. 5). The innermost regions are affected by optical depth effects, the modelled transitions of both OH and OD become optically thick, even before the gas density reaches the critical density ($n_{cr}$).

With respect to the elemental D/H ratio ($\sim 1.5 \times 10^{-5}$, e.g. Vidal-Madjar et al. 1983; Linsky 2003), chemical models predict a considerable enhancement of OD in the diffuse gas (Croswell & Dalgarno 1985), where most of the deuterium is expected to be present in atomic form. A significant amount of OD can then be formed through rapid exchange reactions (Roberts et al. 2002); the reaction D + OH $\rightleftharpoons$ H + OD is exothermic and has a reaction rate coefficient that is higher at low temperatures (Margitan et al. 1975; Atahan et al. 2005; Wang et al. 2006), leading to a stronger enhancement in OD towards colder gas. The destruction of OD mainly happens through reactions with C⁺ ions and photodissociation. Our values for the deuterium fractionation in OD, ranging between $4.8 \times 10^{-3}$ and $2.5 \times 10^{-2}$, are higher than those predicted for the diffuse gas by these chemical models, and are consistent with a picture where in the dense gas, with less C⁺ and UV photons the destruction of OD slows down, leading to a stronger deuterium fractionation compared to the diffuse gas.

## 5.2. Decreasing OD abundance towards the more evolved sources

To investigate the OD abundance variation over our sample of massive clumps in various evolutionary stages, we first used our values for the OD column density, $N_{OD}$, from Table 3 towards all sources with a detection in the $^2\Pi_{3/2}$ OD ground-state line. We estimated the H₂ column density from the 870 µm beam averaged flux density from the ATLASGAL point source catalogue (Csengeri et al. 2014), because the APEX beam at this frequency is comparable to the beam size of 20′, corresponding to the OD observations. For the dust temperature ($T_{dust}$), we used the values estimated by König et al. (2017) based on fitting the spectral energy distribution (SED) of the dust emission. The $T_{dust}$ values are found to be between 22.1 and 34.5 K, the temperatures around 30 K are found towards sources hosting embedded UCHII regions. Then, we computed the observationally estimated OD abundance by $X(OD) = N_{OD}/N_{H_2}$.

We complemented the observationally derived $X(OD)$ estimates (where available) with the results of the radiative transfer modelling (Sect. 4), and hence we can discuss the abundance variations for a sample of nine sources. We find that the highest OD abundances are near the youngest sources of the sample, except for the non-detection towards G34.41+0.2. In particular, G23.21–0.3 and G35.20–0.7 show about an order of magnitude larger abundances than the rest of the sample. Figure 9 shows the distribution of $X(OD)$ as a function of $L_{bol}/M$ that we use as a proxy for the evolutionary stage of the clump (see, however, the work of Ma et al. 2013 about $L_{bol}/M$ reflecting the mass of the central object). Pearson's correlation coefficient ($\rho$) is -0.56, suggesting a trend of decreasing $X(OD)$ over $L_{bol}/M$. Due to the limited number of sources observed, this correlation is, however, not statistically significant (false alarm probability 16%).

We find that the classification of the sources based on the diagnostic of $L_{bol}/M$ does not strictly reflect the evolutionary stage of deeply embedded young stellar objects whose envelopes may just be too compact to impact the $L_{bol}$ on larger scales. An interesting case is G351.58–0.4, which has a low $L_{bol}/M$, and hence it is expected to be among the younger sources. However, as it shows mid-infrared emission, it has been classified
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Observed OD abundance ($X(OD)$) versus $L_{bol}/M$ ratio for the entire sample. The colors show the different source types. The grey arrow marks the upper limit estimated for the non-detection of the quiescent source, G34.41. The Pearson correlation coefficient ($\rho$) is labelled in the figure.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Observed OD and HDO line area measured on the line-to-continuum ratio. The colours show the different source types.</figcaption></figure>

as an infrared bright source in König et al. (2017). Based on radio observations, Culverhouse et al. (2011) find a spectral index consistent with the criteria for UCHII regions from Wood & Churchwell (1989), suggesting that it hosts already formed high-mass star(s). Nevertheless, consistently with its low $L_{bol}/M$, we find a relatively high $X(OD)$, suggesting that the bulk of the gas is chemically young and has not yet been strongly impacted by the star formation process.

Our findings suggest that the OD abundance in the cold dense gas is higher at the onset of the collapse and decreases with time. This is consistent with chemical models predicting that both a higher temperature and stronger radiation field lead to a more efficient destruction of OD. Due to the radiative feedback of newly formed high-mass stars, such conditions occur towards the more advanced evolutionary stages of high-mass star formation. Furthermore, through reactions with O, CO, and C, the OD molecule is destructed in the dense gas; however, these are rather slow processes (Croswell & Dalgarno 1985), suggesting that the lifetime of OD in the cold dense gas is long. The example of G351.58-0.4 suggests that the presence of deeply embedded UCHII regions does not have a strong impact on the OD abundance in the outer layers of the clump, instead the overall $L_{bol}/M$ ratio seems more important.

### 5.3. OD and HDO: photodissociation and constraints on the water chemistry

The formation and destruction of OH (and OD) are tightly linked to the chemical reaction network of water. Since data for the ground-state HDO line is available for roughly half of the sample, for the more evolved sources, we typically compare and correlate here the observed OD and HDO columns in the cold envelope. From the line shapes shown in Fig. 3, we already notice a resemblance suggesting that the absorptions originate from similar conditions within the envelope. We used the measured area for the $^2\Pi_{3/2}$ OD ($J = 5/2-3/2$) line from the line-to-continuum ratio (which is directly proportional to the molecular column density in the case of optically thin emission), and perform the same analysis with the HDO ground-state line. Similarly, as in the case for OD, we fitted two velocity components to the HDO line for W49N and W51e2, and we separated the contribution from the outflow and the $70 \text{ km s}^{-1}$ components, respectively. Figure 10 shows the comparison between the OD and HDO absorption from the ground-state transitions. Our results suggest that the OD and HDO column density is correlated; however, we have to note that this result is dominated by the source, G351.58-0.4, which has by far the largest amount of OD and HDO columns. Unfortunately, due to the lack of HDO observations towards a larger number of sources considered young, we cannot conclude whether such a trend holds. However, a tentative correlation between the HDO and OD column would suggest that OD could be a proxy for HDO, and hence H$_2$O, in the outer envelope.

A correlation between OD and HDO might be expected since the formation pathways are tightly linked between OH and H$_2$O. In the cold gas, OH forms together with water as a result of dissociative recombination from H$_3$O$^+$ (Buhr et al. 2010; van Dishoeck et al. 2014). At low temperatures, however, the high water abundance in ices favours a significant contribution from grain surface processes as well (Tielens & Hagen 1982; Roberts & Herbst 2002; Cuppen et al. 2010). On the other hand, at higher temperatures, OH could form through direct reactions between H$_2$ and O if the energy barrier could be overcome (e.g. Elitzur & Watson 1978; Wagner & Graff 1987; van der Tak et al. 2006). A subsequent reaction with H$_2$ then forms followed by H$_2$O. This reaction would essentially convert all OH to H$_2$O, however, the backward reaction proceeds when a strong UV field is present or a high abundance of atomic hydrogen is prevalent. The destruction of OH and H$_2$O is further differentiated by different reaction channels. While H$_2$O is destroyed by an ion-molecule reaction with C$^+$, OH is destroyed by a neutral-neutral reaction with N (Wakelam et al. 2010). More importantly, the equilibrium of the fundamental production and destruction channel of OH + H$_2$ ⇌ H$_2$O + H depends on the local conditions, such as density and temperature. Altogether, this may lead to differences in the OH and H$_2$O abundances throughout the inner envelope. However, when considering the main production channel for H$_2$O, the same reaction for OD + H$_2$ ⇌ HDO + H is likely to be slower (e.g. Thi et al. 2010). Therefore the OD over OH ratio is expected to change once the chemical reaction pathways change in the inner regions. This is expected to lead to a general decrease of HDO/H$_2$O in the hot inner regions, while the OD/OH abundance is expected to increase, and these differences should become more pronounced with time (Coutens et al. 2014).

Based on our radiative transfer modelling of the OH and OD lines, we constrained the deuterium fractionation for OH in the outer envelope. The best estimates are obtained for G34.26+0.15 and G351.58-0.4, which are $4.8 \times 10^{-3}$ and
$2.5 \times 10^{-2}$, respectively. We obtain an upper limit of the order of $1.5 \times 10^{-2}$ for G35.20–0.7 and a lower limit of $2 \times 10^{-2}$ for G327.29–0.6, which are all sources with low $L_{bol}/M$. Although we do not have a large sample of sources with deuterium fractionation for both OH and H$_2$O, the values collected in Table 5 show that the OD/OH and HDO/H$_2$O values are similar within an order of magnitude. The source G351.58–0.4 shows, however, somewhat larger deuterium fractionation for OD compared to HDO.

While the inner regions are expected to be greatly impacted by the photodissociation of H$_2$O, we cannot put strong constraints on the deuterium fractionation of OD in the hot inner regions here. Our observations exclude, however, a stronger fractionation at regions of a few thousands of au for G34.26+0.15. Observations of higher *J* transitions of OD could put more constraints on the OD/OH fractionation in the inner region, although large optical depths may prohibit these observations.

## 6. Conclusions

We report the detection of the $^2\Pi_{3/2}$ OD ($J = 5/2-3/2$) line at 1.3 THz with the SOFIA/GREAT instrument towards a sample of Galactic high-mass star forming regions. We estimate the OD column densities in a range of $N(\text{OD}) = 4-70 \times 10^{12} \text{ cm}^{-2}$ for these sources, where the highest OD column densities are seen towards the youngest clumps of the sample.

We performed a detailed radiative transfer modelling for one of the best known source of the sample, G34.26+0.15. Our results show that the $^2\Pi_{3/2}$ OD ($J = 5/2-3/2$) absorption originates from the dense and cold outer layers of the envelope. We put constraints on the deuterium fractionation within the inner regions, where, due to sublimation and higher temperatures, gas-phase chemical formation pathways are expected to dominate the chemistry. Although the limits obtained from these models are not stringent, they allow us to exclude a strong enhancement of deuterium fractionation towards the inner envelope.

We estimate the OD molecular abundance in the outer envelope for nine sources of the sample. A correlation with the $L_{bol}/M$ ratio - which is used as a proxy for evolutionary stage - shows a tentative slow decrease in the OD abundance with time. This is consistent with chemical models of the cool, dense gas phase, where the main formation of OD proceeds through rapid exchange reactions leading to a strong enhancement of deuterium fractionation. The OD molecule is then destroyed by slow reactions taking place in the outer envelope. Our results suggest that radiation from the deeply embedded heating sources at this stage does not strongly impact the conditions in the envelope.

We find a deuterium fractionation for OD of ~0.5% towards G34.26+0.15 in the outer envelope, which is consistent with estimates from the HDO/H$_2$O ratio. We find a weak correlation between the OD and HDO column in a sample of 6 sources for which ancillary data were available. This suggests that OD could be a proxy for HDO and thus H$_2$O in these sources.

**Acknowledgements.** We thank the referee for the careful reading of the manuscript and constructive suggestions that improved the paper. Based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy. SOFIA Science Mission Operations are conducted jointly by the Universities Space Research Association, Inc., under NASA contract NASB-97001, and the Deutsches SOFIA Institut under DLR contract 50 OK 0901. T.Cs. acknowledges financial support from the French State in the framework of the IdEx Université de Bordeaux Investments for the future Program. The development and operation of GREAT is financed by resources from the participating institutes, by the Deutsche Forschungsgemeinschaft (DFG) within the grant for the Collaborative Research Center 956 as well as by the Federal Ministry of

Economics and Energy (BMWI) via the German Space Agency (DLR) under Grants 50 OK 1102, 50 OK 1103 and 50 OK 1104.

## References

Atahan, S., Alexander, M. H., & Rackham, E. J. 2005, J. Chem. Phys., 123, 204306

Betz, A. L., & Boreiko, R. T. 1989, ApJ, 346, L101

Buhr, H., Stützel, J., Mendes, M. B., et al. 2010, Phys. Rev. Lett., 105, 103202

Ceccarelli, C., Caux, E., White, G. J., et al. 1998, A&A, 331, 372

Cesaroni, R., & Walmsley, C. M. 1991, A&A, 241, 537

Charnley, S. B. 1997, ApJ, 481, 396

Comito, C., Schilke, P., Gerin, M., et al. 2003, A&A, 402, 635

Coutens, A., Vastel, C., Hincelin, U., et al. 2014, MNRAS, 445, 1299

Croswell, K., & Dalgarno, A. 1985, ApJ, 289, 618

Csengeri, T., Menten, K. M., Wyrowski, F., et al. 2012, A&A, 542, L8

Csengeri, T., Urquhart, J. S., Schuller, F., et al. 2014, A&A, 565, A75

Csengeri, T., Leurini, S., Wyrowski, F., et al. 2016, A&A, 586, A149

Culverhouse, T., Ade, P., Bock, J., et al. 2011, ApJS, 195, 8

Cuppen, H. M., Ioppolo, S., Romanzin, C., & Linnartz, H. 2010, Phys. Chem. Chem. Phys., 12, 12077

Dagdigian, P. J. 2021, MNRAS, 505, 1987

Dalgarno, A. 1994, Adv. Atm. Mol. Opt. Phys., 32, 57

Elitzur, M., & Watson, W. D. 1978, ApJ, 222, L141

Flagey, N., Goldsmith, P. F., Lis, D. C., et al. 2013, ApJ, 762, 11

Fraser, H. J., Collings, M. P., McCoustra, M. R. S., & Williams, D. A. 2001, MNRAS, 327, 1165

Gerin, M., de Luca, M., Black, J., et al. 2010, A&A, 518, L110

Giannetti, A., Wyrowski, F., Brand, J., et al. 2014, A&A, 570, A65

Giannini, T., Nisini, B., & Lorenzetti, D. 2001, ApJ, 555, 40

Goldreich, P., & Kwan, J. 1974, ApJ, 189, 441

Guan, X., Stutzki, J., Graf, U. U., et al. 2012, A&A, 542, L4

Güsten, R., Nyman, L. Å., Schilke, P., et al. 2006, A&A, 454, L13

Harada, N., Herbst, E., & Wakelam, V. 2010, ApJ, 721, 1570

Hatchell, J., Thompson, M. A., Millar, T. J., & MacDonald, G. H. 1998, A&AS, 133, 29

Heyminck, S., Graf, U. U., Güsten, R., et al. 2012, A&A, 542, L1

Higgins, R., Kabanovic, S., Pabst, C., et al. 2021, A&A, 652, A77

Hogerheijde, M. R., & van der Tak, F. F. S. 2000, A&A, 362, 697

Immer, K., Galván-Madrid, R., König, C., Liu, H. B., & Menten, K. M. 2014, A&A, 572, A63

Jensen, M. J., Bilodeau, R. C., Safvan, C. P., et al. 2000, ApJ, 543, 764

Klos, J., Ma, Q., Dagdigian, P. J., et al. 2017, MNRAS, 471, 4249

König, C., Urquhart, J. S., Csengeri, T., et al. 2017, A&A, 599, A139

Kulczak-Jastrzębska, M. 2017, ApJ, 835, 121

Kurayama, T., Nakagawa, A., Sawada-Satoh, S., et al. 2011, PASJ, 63, 513

Lesurf, J. C. G. 1981, Infrared Phys., 21, 383

Leurini, S., Wyrowski, F., Wiesemeyer, H., et al. 2015, A&A, 584, A70

Linsky, J. L. 2003, Space Sci. Rev., 106, 49

Liu, F.-C. 2017, PhD thesis, Rheinische Friedrich-Wilhelms-Universität Bonn, Germany

Liu, F.-C., Parise, B., Wyrowski, F., Zhang, Q., & Güsten, R. 2013, A&A, 550, A37

Ma, B., Tan, J. C., & Barnes, P. J. 2013, ApJ, 779, 79

Margitan, J. J., Kaufman, F., & Anderson, J. G. 1975, Chem. Phys. Lett., 34, 485

Melnick, G. J., Genzel, R., & Lugten, J. B. 1987, ApJ, 321, 530

Mookerjea, B., Vastel, C., Hassel, G. E., et al. 2014, A&A, 566, A61

Neau, A., Al Khalili, A., Rosén, S., et al. 2000, J. Chem. Phys., 113, 1762

Neufeld, D. A., Sonnentrucker, P., Phillips, T. G., et al. 2010, A&A, 518, L108

Nisini, B., Giannini, T., & Lorenzetti, D. 2002, ApJ, 574, 246

Parise, B., Du, F., Liu, F. C., et al. 2012, A&A, 542, L5

Peebles, P. J. E. 1966, ApJ, 146, 542

Peng, T. C., Wyrowski, F., van der Tak, F. F. S., Menten, K. M., & Walmsley, C. M. 2010, A&A, 520, A84

Persson, C. M., Black, J. H., Cernicharo, J., et al. 2010, A&A, 521, L45

Pickett, H. M., Poynter, R. L., Cohen, E. A., et al. 1998, J. Quant. Spectr. Rad. Transf., 60, 883

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C. The Art of Scientific Computing (Cambridge: Cambridge University Press)

Purcell, C. R., Hoare, M. G., Cotton, W. D., et al. 2013, ApJS, 205, 1

Risacher, C., Güsten, R., Stutzki, J., et al. 2016, A&A, 595, A34

Roberts, H., & Herbst, E. 2002, A&A, 395, 233

Roberts, H., & Millar, T. J. 2007, A&A, 471, 849

Roberts, H., Herbst, E., & Millar, T. J. 2002, MNRAS, 336, 283

Roelfsema, P. R., Helmich, F. P., Teyssier, D., et al. 2012, A&A, 537, A17

Rolffs, R., Schilke, P., Wyrowski, F., et al. 2011a, A&A, 527, A68
Rolffs, R., Schilke, P., Zhang, Q., & Zapata, L. 2011b, A&A, 536, A33

Sánchez-Monge, Á., Cesaroni, R., Beltrán, M. T., et al. 2013, A&A, 552, L10

Schuller, F., Menten, K. M., Contreras, Y., et al. 2009, A&A, 504, 415

Storey, J. W. V., Watson, D. M., & Townes, C. H. 1981, ApJ, 244, L27

Thi, W. F., Woitke, P., & Kamp, I. 2010, MNRAS, 407, 232

Tielens, A. G. G. M., & Hagen, W. 1982, A&A, 114, 245

van der Tak, F. F. S., & Menten, K. M. 2005, A&A, 437, 947

van der Tak, F. F. S., Walmsley, C. M., Herpin, F., & Ceccarelli, C. 2006, A&A, 447, 1011

van Dishoeck, E. F., Herbst, E., & Neufeld, D. A. 2013, Chem. Rev., 113, 9043

van Dishoeck, E. F., Bergin, E. A., Lis, D. C., & Lunine, J. I. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 835

van Dishoeck, E. F., Kristensen, L. E., Mottram, J. C., et al. 2021, A&A, 648, A24

Vidal-Madjar, A., Laurent, C., Gry, C., et al. 1983, A&A, 120, 58

Wagner, A. F., & Graff, M. M. 1987, ApJ, 317, 423

Wakelam, V., Herbst, E., Le Bourlot, J., et al. 2010, A&A, 517, A21

Walsh, A. J., Burton, M. G., Hyland, A. R., & Robinson, G. 1998, MNRAS, 301, 640

Wampfler, S. F., Bruderer, S., Kristensen, L. E., et al. 2011, A&A, 531, L16

Wampfler, S. F., Bruderer, S., Karska, A., et al. 2013, A&A, 552, A56

Wang, W., Santos, E., & Brandão, J. 2006, J. Chem. Phys., 124, 074305

Weaver, H. F., & Williams, D. R. W. 1964, Nature, 201, 380

Weinreb, S., Barrett, A. H., Meeks, M. L., & Henry, J. C. 1963, Nature, 200, 829

Wienen, M., Wyrowski, F., Schuller, F., et al. 2012, A&A, 544, A146

Wiesemeyer, H., Güsten, R., Heyminck, S., et al. 2012, A&A, 542, L7

Wiesemeyer, H., Güsten, R., Heyminck, S., et al. 2016, A&A, 585, A76

Wilson, T. L. 1999, Rep. Prog. Phys., 62, 143

Wolfire, M. G., & Cassinelli, J. P. 1986, ApJ, 310, 207

Wood, D. O. S., & Churchwell, E. 1989, ApJS, 69, 831

Wyrowski, F., Güsten, R., Menten, K. M., Wiesemeyer, H., & Klein, B. 2012, A&A, 542, L15

Wyrowski, F., Güsten, R., Menten, K. M., et al. 2016, A&A, 585, A149

Yang, Y.-L., Green, J. D., Evans, II, N. J., et al. 2018, ApJ, 860, 174

Young, E. T., Becklin, E. E., Marcum, P. M., et al. 2012, ApJ, 749, L17
## Appendix A: Model results for G351.58–0.4

In Fig. A.1, we show the results of our RATRAN modelling for G351.58–0.4. As for the other sources, these models assume a constant abundance profile for all species. While the fit for the OD $^2\Pi_{3/2}$ $J=5/2-3/2$ transition at 1390.61 GHz is satisfactory, the OH $^2\Pi_{1/2}$ $J=3/2-1/2$ transition is not well reproduced by this model. As suggested by the NH$_3$ ($J=3_2-2_2$) line profile from Wyrowski et al. (2016), the blueshifted component of a broader absorption component is apparent, which also seems to impact the excited-state OH line profiles. Furthermore, our assumption on the constant OH abundance profile may also be inadequate for this source.

## Appendix B: Comparison with the NH$_3$ ($J=3_2-2_2$) 1.8 THz absorption

We compare the line profiles of the OD $^2\Pi_{3/2}$ $J=5/2-3/2$ line with the NH$_3$ ($J=3_2-2_2$) line in Fig. B.1. For better visibility, the NH$_3$ lines are scaled. We notice that the line profiles are remarkably similar, showing very similar rest velocities. Towards the sources G327.29–0.6 and G351.58–0.4, absorption from a blueshifted broad velocity component, likely due to molecular outflows, is visible in the NH$_3$ line, while it is absent in the ground-state OD line.

## Appendix C: Removal of instrumental artefacts from the $^{18}$OH $^2\Pi_{3/2}$ ground-state line

In supra-THz heterodyne spectroscopy, local oscillator (LO) power is precious, and transmission losses in the injection of the generated signal are to be kept at the minimum. The $^{18}$OH $^2\Pi_{3/2}$ ground-state line was obtained in 2013, when the LO signals of GREAT (Heyminck et al. 2012) were injected by means of Martin-Pupplet diplexers (e.g. Lesurf 1981). Depending on the tuning, standing waves may arise, causing sinusoidal baselines that are further modulated by the limited frequency passband of the diplexer. Moreover, beyond standing wave features, gain fluctuations can originate in the mixing or subsequent amplification stages. If they are faster than the occurrence of calibration scans, they further deteriorate the baseline quality, especially in observations towards far-infrared, strong continuum sources. These effects together impact the tuning of the Ma band to the $^{18}$OH $^2\Pi_{3/2}$ line.

By their nature, these standing waves cannot be removed by a simple Fourier filter, unlike their equivalents occasionally arising between the mixer and the telescope's secondary mirror. The fitting of high-order polynomial baselines (~10) was found to introduce spurious line profiles at the systemic velocity and is, likewise, to be avoided. Adjusting sinusoidal waveforms, although numerically somewhat more stable, lacks objectiveness as well. In contrast, differential total power spectra, extracted from the off-source phases of subsequent chop-and-nod sub-scans and typically separated by ~1 min, were observed to exhibit similar baseline distortions while being void of the astronomical signal (Higgins et al. 2021, further references therein). In the stable part of the spectral bandpass the aforementioned artefacts are linear in nature, which suggests the necessity to use a linear combination of such differential off-source total power spectra (hereafter $a_j$, $j=1,2,\dots,n$) to fit the residual baseline apparent in the averaged on-off spectrum, denoted **b**. All spectra are represented in an *m* dimensional vector space, where *m* is the number of spectral channels considered. For this purpose,

we minimised the quantity:

$$ \chi^2 = |A \cdot \mathbf{p} - \mathbf{b}|^2, \qquad (C.1) $$

where the columns of the $m \times n$ matrix A are the spectra $a_j$, and where the coefficients of the linear combination minimizing $\chi^2$ are denoted as a vector $\mathbf{p} = (p_1, p_2, \dots, p_n)^T$. Since $m > n$, the problem to be solved is over-determined; yet, a direct minimisation using normal equations is often numerically unstable, because ambiguous combinations of coefficients $\mathbf{p}$ can occur, causing the Jacobian to become singular. However, the solution vector $\mathbf{p}$ can be obtained from the easily invertible singular value decomposition (SVD) of matrix A, providing a more stable method to minimise $\chi^2$ (e.g. Numerical Recipes, Press et al. 1992). Our baseline fit can then be constructed as $\mathbf{b} = \sum_{j=1}^n p_j a_j$.

In practice, we do not directly use the spectra $a_j$, we used smoothed fitting functions instead; this is to avoid excess noise otherwise introduced by the baseline removal. Third order wavelet transforms proved to provide baseline fits with the smallest $\chi^2$. Other similarly smoothed variants may yield satisfactory results alike, while higher orders bear the risk of fitting the spectral line as a baseline feature even though the spectra $a_j$ do not contain the astronomical signal. As a safeguard against this, we use a weighted fit, modifying Eq. C.1:

$$ \chi^2 = \left|\left(\mathrm{diag}(w_j)\right) \cdot (A \cdot \mathbf{p} - \mathbf{b})\right|^2, \quad j = 1, \dots, m, \qquad (C.2) $$

with weights

$$ w_j = 1 - \exp\left(-4 \ln 2 \cdot \left(\frac{v_j - v_0}{\Delta_v}\right)^2\right), \qquad (C.3) $$

centred at $v_0 = 60$ km s$^{-1}$. The exact choice of the width $\Delta_v$ is not critical as long as the putative $^{18}$OH absorption line is covered, yet leaving enough spectral baseline on either side (i.e., $\Delta_v = 10$ to $15$ km s$^{-1}$). For the computations we use *Numerical Recipes* routines svdcmp and svbksb to perform the SVD and to calculate $\mathbf{p}$, respectively.

The application of the SVD technique to the baseline removal from the $^{18}$OH $^2\Pi_{3/2}$ line is demonstrated in Fig. C.1. After iterative Hanning smoothing to 0.6 km s$^{-1}$ channel separation, the resulting spectrum features a sensitivity of $\sigma_{rms}(T_A^*) = 207$ mK prior to transmission correction, which equals, within errors, the radiometric noise pertaining to the measured single-sideband system temperature in the relevant part of the bandpass, 6070 K. The obtained spectrum therefore qualifies to define a lower limit to the $^{16}$OH/$^{18}$OH ratio (Sect. 4 and Fig. 6).

## Appendix D: Comparison of RATRAN models with different collisional rate coefficients

In Sect. 4, we used RATRAN to model the OD absorption, and adopted the collisional rate coefficients of OH (Kłos et al. 2017) also for OD. However, recently, collisional rate coefficients for OD with H$_2$ became available (Dagdigian 2021) that also include the hfs splitting. To test the validity of our initial approach, we compared our best fit model for the OD 1391 GHz line (Sect. 4.2) for G34.26+0.15 both using the OH coefficients from Kłos et al. (2017) adapted for OD, and the collisional rate coefficients of OD including the hfs structure from Dagdigian (2021). We show this comparison in Fig. D.1 and see that the difference between the models (red lines) is marginal.
<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Same as Fig. 6, but for G351.58−0.4. The observed spectra are shown in black, red line shows the fitted RATRAN model. Dotted line marks the $v_{lsr}$ of the source. The panel for the $^2\Pi_{3/2}$ J=5/2-3/2 transition at 2514 GHz shows the line-to-continuum ratio ($T_1/T_c$), all other panels show the line temperature in Kelvins.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. B.1.** Comparison between the OD $^2\Pi_{3/2}$ ($J = 5/2-3/2$) transition (black) and the NH$_3$ ($J = 3_{2+}-2_{2-}$) transition (red).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. D.1.** Comparison of our RATRAN models for the OD 1391 GHz line using the collisional rate coefficients for OD from Dagdigian (2021) (left) and the collisional rate coefficients from OH adapted for OD from Klos et al. (2017) (right).</figcaption></figure>

<caption>**Fig. C.1.** Spectrum of $^{18}$OH$^2\Pi_{3/2}$ before (top) and after (bottom) removal of the SVD derived baseline (red line, top) and transmission correction. The relatively narrow kernel of a telluric O$_3$ line is deliberately left uncorrected. The blue line represents the weight profile used in the $\chi^2$ minimisation.</caption>